Finite-Temperature Screw Dislocation Core Structures and Dynamics in α -Titanium

A multiscale approach based on molecular dynamics (MD) and kinetic Monte Carlo (kMC) meth-ods is developed to simulate the dynamics of an ⟨ a ⟩ screw dislocation in α -Ti. The free energy barriers for the core dissociation transitions and Peierls barriers for dislocation glide as a function of temperature are extracted from the MD simulations (based on Machine Learning interatomic potentials and optimization); these form the input to kMC simulations. Random walk dislocation trajectories from kMC agree well with those predicted by MD. On some planes, dislocations move via a locking-unlocking mechanism. Surprisingly, some dislocations glide in directions that are not parallel with the core dissociation direction. The MD/kMC multiscale method proposed is applicable to dislocation motion in simple and complex materials (not only screw dislocations in Ti) as a function of temperature and stress state.


I. INTRODUCTION
The major plastic deformation mechanism in crystalline metals is dislocation glide.The motion of dislocations is (largely) controlled by elastic stresses and intrinsic dislocation properties.While these stresses are easily analyzed in terms of continuum elasticity, dislocation motion often occurs preferentially on planes other than those with the largest resolved shear stress.This is associated with differences in glide resistance/lattice friction between different slip planes.This depends on the relative ease with which the dislocation core moves.The ease of glide, in turn, is sensitive to the dislocation core structure [1][2][3][4][5][6].The dislocation core structure is, in many cases, temperature-dependent.While knowledge of the core structure is essential, a quantitative link between core structure and dislocation dynamics is often elusive.Here, we develop a multiscale approach to predict screw dislocation dynamics in α-Ti.
Atomistic simulations are commonly employed to determine dislocation core structure, in part because the direct experimental determination of the structure is demanding [7].Since the atomic structure of the material is highly distorted with respect to that in a perfect crystal, quantum mechanical accuracy is often required to predict core structures; often achieved using density functional theory (DFT) calculations [8][9][10][11].Transition state theory-based methods (such as the nudged elastic band NEB method) are often employed to discern the minimum energy path of a dislocation core as it traverses the slip plane [10][11][12][13].While DFT methods are usually limited to ground-state (0 K) structures, finite-temperature ab initio molecular dynamics are possible but too computationally costly for widespread use.Therefore, most finite-temperature dislocation core structure determination is based upon semi-empirical interatomic potential methods; e.g., Poschmann et al. [14] studied the core structure of an ⟨a⟩ = a⟨1010⟩ screw dislocation in α-Ti at finite temperature using a modified embedded atom method (MEAM) potential [15].Unfortunately, this MEAM potential fails to accurately reproduce all of the relevant 0 K core structures and energies predicted by DFT [16].Bond order potentials (BOPs) [17,18] were proposed to retain the quantum nature of atomic interactions in transition metals in a more cost-effective manner than DFT.However, BOPs are both computationally costly and not easily implemented in molecular dynamics (MD) simulations [19].The recently developed Deep Potentials (DPs) [20,21] (a class of neural network potential) yield DFT accuracy with near empirical potential computational efficiency.Here, we employ the DP method to predict the structure of screw dislocation cores at finite-temperature in α-Ti.
The link between dislocation core structure and dislocation dynamics is often related to the assumption that the dislocation glide direction is consistent with the dislocation core dissociation direction.While such an assumption may be valid in some simple cases, its validity is far from assured in the case of more complex (non-cubic) materials, such as hexagonal close packed (HCP) metals.Several models have been proposed to simulate dislocation dynamics at the mesoscale, such as discrete dislocation dynamics (DDD) [22][23][24] and kinetic Monte Carlo (kMC) [23,[25][26][27][28] methods.Extant mesoscale models do not explicitly incorporate the effects associated with the dislocation core structure.Here, we develop a mesoscale dislocation dynamics model that incorporates an explicit description of the atomic-scale character of dislocation core structure.
While recent simulations (e.g., see [14]) focus on long dislocation lines, here we focus on the intrinsic dislocation properties associated with the dislocation core structure.The admittedly important roles played by dislocation kinks in the motion of long dislocations are omitted here in order to provide a thorough examination of core effects without the complicated features of kink dynamics (which vary dramatically with, e.g., local dislocation curvature, junctions and interactions with other dislocations).The effects of core structure of short ⟨a⟩ screw dislocation segments in α-Ti based on the DP for Ti as a function of temperature are investigated.Experimental observations [29] show that the edge dislocations are highly mobile and the yield strength of HCP Ti is governed by screw dislocation lattice friction.We report the results of MD simulations (based on machine learning potentials of quantum mechanical accuracy) of screw dislocation core structures, transitions between different core structures, statistical analysis of dislocation motion, and the determination of the kinetic parameters (i.e., free energy barriers associated with migration, core structure transitions, ...) describing screw dislocation motion.We then perform kinetic Monte Carlo simulations of screw dislocation motion in Ti incorporating these quantum mechanically accurate MD simulation parameters.We examine the effect of both temperature and loading direction on dislocation core transitions and dislocation mobility.The results provide the basis for understanding non-Arrhenius screw dislocation mobility in metals with complex crystal structures.

II. ATOMISTIC DISLOCATION STRUCTURE AND DYNAMICS
A. Dislocation Core Structure We performed MD simulations to determine ⟨a⟩ screw dislocation core structure in α-Ti; the simulation geometry is shown in Fig. 1a (see Methods).By examining all the MD configurations from 300-900 K, we identified five distinct core structures for the ⟨a⟩ screw dislocation in HCP Ti. Figure .1c-g show these five core structures with a differential displacement map [30] (the black arrows) and the Nye tensor component α zz [31,32] (the contour).The α zz map describes the screw component of the Burgers vector density.We find that the distribution of α zz is highly delocalized in the form of a dipole, indicating that the core structure dissociates on a plane that includes [1 210].For convenience, we denote the pyramidal plane by "π", the prismatic plane by "P" and the basal plane by "B".We find that two of the dislocation cores dissociate along the π plane; i.e., the "π core" (Fig. 1c) and "π ′ core" (Fig. 1d).The dissociation plane of π core is close to ( 3031) while that of π ′ core is close to ( 1011), as shown in Fig. 1b.Two cores are dissociated along the P plane; i.e., the "P core" (Fig. 1e) and "P ′ core" (Fig. 1f).The α zz (x, y) map for the P core possesses inversion symmetry roughly about the point (0, 0) while that for the P ′ core possesses mirror symmetry roughly about the y = 0 line.We also identify a dissociated core along the B plane; i.e., the "B core" (Fig. 1g).The π, π ′ , P and P ′ cores were found at all temperatures in our simulations.These three are consistent with the 0 K core structures predicted by DFT calculations [10].The B core was observed only at high temperature, T ≳ 400 K.
Since the π, π ′ , P and P ′ cores are stable at 0 K, we can obtain their 0 K equilibrium structures by direct energy minimization based on different initial configurations (with the dislocation core centered at different positions).The π, π ′ , P and P ′ core energies are E π = 544.8± 0.43 meV/ Å, E π ′ = 561.5 ± 0.52 meV/ Å and E P = E P ′ = 547.4± 0.34 meV/ Å, respectively (see Supplementary Information, SI, for details).The energy differences are around E P − E π = 2.6 meV/ Å and E π ′ − E π = 16.7 meV/ Å; for comparison, the DFT re- sults [10] are 5.7 meV/ Å and 11 meV/ Å, respectively.The dislocation core structures and energies at 0 K obtained from DP are reasonably consistent with DFT results.Since the π ′ core energy is much higher than the energies of the other cores and a nudged-elastic-band calculation [10] shows that the π ′ core energy is almost as high as the barrier for π core glide, π ′ core is not important for the thermodynamic and kinetic properties; hence, we ignore the π ′ core below.Since the P and P ′ core energies are nearly equal, we do not distinguish the P and P ′ cores below.

B. Dislocation Core Dynamics
We focus on two aspects of dislocation core dynamics: core motion and core structure transitions.The former is analyzed based on the trajectory of the core position (Fig. 2a), while the latter requires recognition of instantaneous core structure.We automated the determination of the core position and the core dissociation direction based on the α zz map.For convenience, we denote the α zz -weighted average of any quantity A as Ā = A(x, y)α zz (x, y)dxdy α zz (x, y)dxdy , Then, the core position is r core ≡ (x, ȳ); the core positions determined in this manner are indicated by the blue "+" in Figs.1c-g. Figure 2a shows the core trajectory at 500 K. From the trajectory, we see that the dislocation random walk is anisotropic.The pattern is elongated along the y-axis, suggesting that the dislocation glides on the P plane most frequently.
The glide direction as a function of time ϕ(t) is where ∆t is the time step.
Figure 2b shows the temporal evolution of the dissociation direction θ at 500 K and Fig. 2c shows the distribution of the core dissociation direction θ at different temperatures.The three sharp peaks at θ≃75 • , 90 • and 105 • correspond to the π 1 , P and π 2 cores (π 1 and π 2 are both π cores).The peaks for the π cores shift towards θ = 90 • as temperature increases because the c/a ratio increases with temperature [33].The peak for the π ′ core, if it existed, would be at ∼ 61 • and ∼ 110 • ; however, no peaks exist there, implying that the π ′ core may be ignored.At a high temperature (e.g., 900 K), there are shallow and broad humps at θ = 0 or 180 • , indicating the existence of the B core at high temperature.The peaks, signaling different cores, broaden with increasing temperature.We set criteria to distinguish the core structures based on the distribution in Fig. 2c.We define the width of the orientation window for the P core as the minima between the P and π cores; i.e., ∼ 90 ± 9 • (the exact position of the minima varies with temperature); see the two dashed lines near θ = 90 • in Fig. 2b or c.The boundary between the π and B cores is not well defined.Since the probability for θ < 50 • or > 130 • is almost zero, any quantity evaluated based on the core distribution is insensitive to the choice of this boundary.In practice, we choose the boundary between the π and B cores at the θ in the middle between 0 and the θ for the π core (i.e., the θ for the second highest peak); see the two dashed lines close to θ = 90 • in Fig. 2b or c.Using this criterion, we can identify the core structure at any time during the MD simulation.We color each point on the trajectory (at 500 K) shown in Fig. 2a according to the core structure at each time.We observe that B cores (green points) are rare.The π core (red points) and the P core (blue points) are largely distributed on alternating P planes; this is consistent with the examination of the π and P core positions in the Fig. 2a inset.The Fig. 2a inset shows that ideally the π core is positioned between a dark gray P plane and a light gray P plane, while the P core and B core are between two light gray or two dark gray P planes.Hence, the π and P cores should be distributed on different P planes.Such consistency validates the core structure recognition method.
It is usually assumed that the dislocation core dissociation direction is the same as the dislocation glide direction; this need not be true.To investigate the correlation between the core dissociation direction (θ) and the glide direction (ϕ), we examine the probability of different glide directions for particular core dissociations (characterized by the dissociation direction θ i (i ∈ {π, P, B}). Figure 2d shows the glide direction distribution for the B core.Most frequently, the B core glides on the B plane (B-glide), corresponding to the peak at |θ − ϕ B | = 0.The B-glide of the B core is also seen in the difference between the α zz maps at a pair of times: ∆α zz ≡ α zz (t 2 ) − α zz (t 1 ).As shown in Fig. 2g, the alternatively distributed negative and positive ∆α zz clouds lying on the B plane is a feature of B-glide of B core.The glide direction distribution for the π core is shown Fig. 2e.This distribution exhibits a peak at |θ − ϕ π | ≈ 15 • (at T = 300 K), where θ π ≈ 75 • .The angle between the π ′ plane (i.e., the ( 1011) plane) and the B plane at 300 K is If the π core glides on the P plane (P-glide), |θ − ϕ π | ≈ 15 • .This sug- gests the peak at ∼15 • in Fig. 2e has contributions from both P-glide and π-glide.Careful examination, however, shows that P-glide is much more frequent than π-glide.
The π-glide and P-glide of π core can be directly verified by the ∆α zz maps in Fig. 2h and i.The alternating negative and positive ∆α zz clouds on the π plane is a feature of π-glide of the π core while the off-line distribution of the negative and positive clouds is a feature of P-glide of the π core.The observation that a π core can glide on the P plane contradicts the assumption that the core glide and the core dissociation directions must be the same.Figure 2f shows the glide direction distribution for the P core.Clearly, the P core only glides on P plane.Again, the P-glide of the P core can be verified through considerations of the ∆α zz clouds distributed along the P plane, as shown in Fig. 2j.

III. MODEL AND PARAMETERIZATION OF DISLOCATION CORE DYNAMICS A. Kinetic Events
The unit kinetic event during the motion of a dislocation core observed in the MD simulation can be abstracted as a transition.We label the states before and after a transition by i and j, respectively, where i, j ∈ {B, π 1 , π 2 , P}.The π 1 and π 2 cores denote, respectively, the π cores with the dissociation directions θ π ≈ 75 • and 105 • (they are symmetry-related).We use κ to label the slip plane, i.e., κ ∈ {B, π 1 , π 2 , P, 0}, where the π 1 and π 2 planes are, respectively, the π ′ planes with the inclination angles ϕ ≈ 61 • and 119 • , and κ = 0 denotes the transition which does not involve the change in core position.We denote the transition from an i to j core involving glide on the κ plane as "i(κ)j".An "i ̸ = j" event represents a transition of the core structure, while an "i = j" event represents glide of a core with no core structural transition.All possible kinetic events observed in the MD simulations (Fig. 2) are shown schematically in Figs.3a-f.These events, denoted i(κ)j, are summarized in Fig. 3g.Some events listed in the 1st, 2nd and 3rd columns are equivalent.The 4th column shows the irreducible events.Note that π denotes the core symmetrically related to π; e.g., if π = π 1 , then π = π 2 .

B. Dislocation Core Free Energy
From the MD results shown in Fig. 2c, we can extract the equilibrium probabilities {P i } that a core is of type i as a function of temperature.The open symbols in Fig. 4a represent the MD data {P i }. π is the most probable core structure at all temperatures,P π (T ) > 0.6.At T < 300 K, the B core is not observed in the MD and at T > 300 K, the B core occurs with very low probability.
In thermal equilibrium at temperature T , the probability of finding the i core is well-described by a Boltzmann distribution: P i ∝ exp (−F i ℓ/k B T ), where F i is the free energy of the i core (energy per length), ℓ is the length of the dislocation line, and k B is the Boltzmann constant.The free energy difference between the i and j cores is related to their probability ratio: TABLE I. Fit parameters for the free energy difference between cores i and j, ∆Fij, as per Eq. ( 5).
3.24 ± 0.266 3.59 ± 0.506 Since the elastic energy is the same for all core structures, ∆F ij represents the core energy difference.The core energy difference ∆F ij at each temperature can be obtained from the MD data {P i } (reported in Fig. 4a) and Eq. ( 4).
Note that although the 10-20 step partial relaxations reduce the potential energy of the system, these energies are not of interest.Rather the partial relaxation simply aids the identification of the inherent core structures.
The important energy differences ∆F ij are determined based on the relative probabilities of different dislocation core structures -which are unaffected by the partial relaxations.The core energy differences, ∆F Pπ and ∆F Bπ , are shown as open symbols in Fig. 4b.We can fit the data ∆F ij with an empirical relationship of the form: where A, B and C are parameters.The rationale for the form of this fitting relation is discussed in the SI.For ∆F Pπ , A = E P − E π = 2.6 meV/ Å, where E P and E π are the energies of the P and π cores at 0 K.The fitted curves are shown as solid lines in Fig. 4b and the parameters are in Table I.

C. Core Dynamics Parameterization
The basic kinetic parameters for dislocation core dynamics are the frequencies of all events as a function of temperature.Unfortunately, it is impractical to deduce these frequencies directly from the MD results.Atomic vibrations are inevitable in MD simulations at finite temperatures; removing these by thermal averaging or quenching in order to unambiguously recognize each core event (defined in Fig. 3) requires artificial criteria.Sampling frequency also makes this impractical: if the sampling frequency is too high, the thermal vibration issue will be severe; if it is too low, we will miss some kinetic events.Here, we sidestep the frequency issue, as explained in this section.In short, we fit the MD data using harmonic transition state theory (HTST) with temperature-independent parameters.The general parameterization steps are S1: (Prediction) Guess a set of frequencies for all kinetic events {ν i(κ)j } at each temperature (see Methods).FIG. 3. Kinetic events for ⟨a⟩ screw dislocation core dynamics.(a)-(f) All possible events starting from a P core, a π core or a B core.The dark gray and light gray circles denote atoms on successive (1 210) planes.The vertical gray solid lines denote a series of ( 2020) planes; the vertical gray dashed lines are offset with respect to the vertical solid lines by √ 3a/4.The horizontal gray solid lines denote a series of (0001) planes (N c planes); the horizontal gray dashed lines are offset with respect to the horizontal solid lines by c/2 (i.e., N c/2 planes).The "+" symbols denote dislocation core positions and their colors are consistent with those in the first column of (g).When two "+" symbols are located at the same site, we make one larger than the other for clarity.(g) The 1 st and 3 rd columns show the starting and ending core structure for one event.The 2 nd nd column shows the core displacements corresponding to the transitions.Some events listed in the 1 st , 2 nd and 3 rd columns are connected by the gray lines indicating that they are equivalent in the sense that they have the same energy landscapes.The 4 th column shows the irreducible events, where i(κ)j denotes the transition from an i to j core by glide on plane κ.
to compute the core mean squared displacement (MSD), the mean squared angular displacement (MSAD) of the core dissociation direction, and the probability of occurrence of each core structure {P i }.
S3: (Optimization) Optimize {ν i(κ)j } such that the MSD, MSAD and {P i } obtained from the kMC are consistent with the MD results at all temperatures (see Methods).
The mean squared displacement MSD and mean squared angular displacement MSAD (S2 and S3 ) are defined as follows.At a particular temperature, the mean squared displacement MSD in the e x -and e y -directions are defined as where ⟨•⟩ denotes the average over t.The circles in Figs.5a and b show the MSDs, ⟨(∆x) 2 ⟩ and ⟨(∆y) 2 ⟩, obtained from the MD simulations under different temperatures.We find that core glide along the e y -axis (P plane) is, in general, faster than glide along the e x -axis (B plane).At all temperatures, each MSD is approximately a linear function of τ .The translational diffusion coefficients in the e x -and e y -directions are One goal of optimizing the frequencies {ν i(κ)j } is to ensure that the values of D T x and D T y from kMC and MD simulations match.The mean squared angular displacement MSAD for core dissociation i is where the subscript "i" denotes that in the average θ(t = 0) = θ i , i.e., the dissociation direction of core i.The circles in Figs.5c, d and e show the MSADs, ⟨(∆θ) 2 ⟩ P , ⟨(∆θ) 2 ⟩ π and ⟨(∆θ) 2 ⟩ B , obtained from the MD simulations at different temperatures.The MSAD, at each temperature, is well fitted by the function: where D R i (i = P, π, B) is the rotational diffusion coefficient about dissociation angle θ i [34].In S3, the frequencies {ν i(κ)j } are trained for the best match between the core probabilities ({P i }), the translational diffusion coefficients (D T x and D T y ) and the rotational diffusion coefficients (D R P , D R π and D R B ) obtained by kMC and the MD results.The solid lines in Figs.5a-c show the MSDs and MSADs obtained from kMC simulations with optimized {ν i(κ)j } for different temperatures.The kMC and MD results are in excellent agreement.
We now turn to the parameter space reduction in S4.The frequencies {ν i(κ)j } are obtained via S1-3 for the MD simulation temperatures.In principle, dislocation core dynamics at other temperatures may also be obtained from MD simulation at such temperatures and repeating S1-3.The computational resources required for these MD simulations and the optimization process limit the applicability of this approach.We resolve this issue based upon a set of additional assumptions.
Consider the schematic energy landscape for a kinetic event i(κ)j in Fig. 6a.Two local minima correspond to the i and j core free energies, F i and F j .The core energy difference is ∆F ij (Eq.( 4)); for a core glide event, ∆F ii = 0.The total free energy barrier for i(κ)j is denoted F b i(κ)j .The intrinsic free energy barrier for i(κ)j is Q i(κ)j and, in principle, Q i(κ)j = Q j(κ)i .The total free energy barrier is commonly approximated as where the factor 1/2 is valid when ∆F ij ≪ Q i(κ)j [26].
While other reasonable proposals for F b i(κ)j are possible, our fitting results, below, suggest that Eq. ( 10) reproduces the MD results.
The frequency of event i(κ)j can be expressed as based on HTST [35,36], where ν 0 i(κ)j is an attempt frequency that includes the effect of barrier recrossing.Equation ( 11) can be rewritten as This shows that η i(κ)j is a function of the inverse temperature T −1 .If {ν 0 i(κ)j } and {Q i(κ)j } are constant with respect to temperature, they can be obtained by linear fitting.However, we find that the fitting quality is improved by allowing the intrinsic free energy barrier to be temperature-dependent and of the form: where Q 0 i(κ)j is the intrinsic barrier at 0 K, T 0 = 1250 K is the HCP-BCC transition temperature for this DP potential [33], and we assign q = 3 to give best fit to all the data.Equation ( 13) is proposed based on the observations that (i) the free energy barrier decrease with increasing temperature and (ii) when the HCP phase becomes unstable/metastable, the transition between the core structures is barrierless.In this way, the kinetic parameters can be fitted to the MD data; i.e., {ν 0 i(κ)j } and {Q 0 i(κ)j }.
Figures 6b and c show the fitting results for Eqs. ( 12) and (13) (symbols/lines are the MD data/fits).As expected, the i(κ)j and j(κ)i data coincide.Among the transition events (Fig. 6b), the transitions between P and π cores are the most frequent and associated with the lowest energy barriers.As expected, direct transitions between P and B cores are rare.Among the glide events (Fig. 6c), the glide on B plane is associated with the lowest barrier.But glide on B plane is rare since the probability of the B core is very low (Fig. 4a).Pyramidal glide (via π(π)π) is much less frequent and asso- Fi and Fj are the free energies of the i and j cores; ∆Fij = Fj − Fi and ∆Fij = 0. F b i(κ)j is the total free energy barrier and Q i(κ)j is the intrinsic barrier for i(κ)j.(b) η i(κ)j (Eq.( 12)) vs. inverse temperature.(c) η i(κ)i (Eq.( 12) for glide events) vs. the inverse temperature.In (b) and (c), symbols and solid lines denote MD data and fitting results (Eqs.( 12) and ( 13)).ciated with a much higher barrier than prismatic glide (via P(P)P).This means that pyramidal glide is rarer than prismatic glide; in qualitative agreement with the DFT 0 K glide barriers [10] (DFT predicts the Peierls barrier for P(P)P glide and π(π)π glide as 11.4 meV Å−1 and 0.4 meV Å−1 , while the values from our simulations are 14.1 meV Å−1 and 1.14 meV Å−1 ).
TABLE II.The kinetic parameters determined by fitting to the MD data for the core transition events and glide events.ν 0 ik is the attempt frequency in Eq. ( 12).Q 0 ik is the intrinsic 0 K energy barrier in Eq. ( 13).

A. Random Walk of a Dislocation Core
The equilibrium core structure probabilities {P i } were computed via kMC and compared with those obtained from MD, as shown in Fig. 4a.At temperatures ≥ 300 K, the kMC results show excellent agreement with MD.P π and P P exhibit a minimum and maximum near 300 K. Beyond this temperature both P π and P B increase, while P P decreases with temperature T .P P and P B are approximately equal near 0 and 900 K.No MD is available below 300 K where it is difficult to obtain valid statistics; here, we only show kMC data from Eq. ( 5).At 0 K, based on the fact that π core is energetically favorable, P π should be 1 and P P and P B should be 0. P B ≈ 0 for T ≤ 300 K, which is consistent with the MD observation that the B core is unstable.
The MSDs and MSADs computed via MD and kMC simulations are compared in Fig. 5a-e.In general, the kMC simulations reproduce the MD results well, except for the MSADs corresponding to core dissociation starting from θ B at low temperatures (Fig. 5e) (this is associated with limited sampling of the B core in the MD simulations).Due to larger accessible timescale in kMC compared with MD, the kMC simulations provide smooth curves at low computational cost (compared with MD).The glide-direction dependent translational diffusion co-efficient can be constructed as see Fig. 5f.D T (ϕ) is elongated in the e y -direction, indicating that a dislocation core moves fast along the P plane and slowly along the B plane; this is consistent with the MD trajectories in Fig. 2a.The dissociationangle-dependent rotational diffusion coefficient can be constructed as where (i, j) = (B, π) or (π, P); see Fig. 5g.D R (θ) measures the rotation rate for cores initially oriented at angle θ. Figure 5g shows that D R (θ) is highly anisotropic at low temperatures and becomes more isotropic as temperature increases.At all temperatures (except for 900 K), D R is maximum at θ = 0 which corresponds to the B core.This is consistent with the fact that B core is not energetically favorable and tends to transform into other cores.D R is a minimum at ∼θ = 75 • ; corresponding to the π core.This indicates that the most stable core is π.

B. Stress-Driven Dislocation Motion
The kMC model is sufficiently flexible to simulate dislocation dynamics under an externally applied stress.We apply the approach developed by Ivanov and Mishin [37].
The major assumption is that the applied stress influences the dislocation glide barrier through the resolved shear stress (RSS), but not the energy barrier for the core transition itself.Hence, this model does not fully capture the non-Schmid effect [38][39][40][41] (see below).
We assume that the free energy landscape for i core glide on the κ plane has the form of where r is a slip distance, Q i(κ)i is the glide barrier (Eq.( 13)) and L i(κ)i is a lattice period in the slip direction on the κ plane.Applying an external stress σ creates Gibbs free energy landscape G i(κ)i (r) = F i(κ)i (r)−f κ r, where f κ is the Peach-Koehler (PK) force: With our sample geometry, Fig. 1a, ξ = e z is the line direction, b = be z is the Burgers vector, s κ = cos ϕ κ e x + sin ϕ κ e y is the slip direction on the κ plane with inclination angle ϕ κ , and τ κ ≡ −σ xz sin ϕ κ + σ yz cos ϕ κ is the RSS on the κ plane.Then, the Gibbs free energy barrier for the forward/backward glide is where τ 0 i(κ)i ≡ πQ i(κ)i /(L i(κ)i b).The frequency of forward/backward glide event i(κ)i, ν ± i(κ)i , is obtained from Eq. ( 11) with . The detailed explanation of Eq. ( 16) can be found in Ref. [37].With this frequency as input, we perform kMC simulations of dislocation motion under different stresses σ (see Methods).
The kMC simulations show a locking-unlocking type of dislocation motion at low temperatures.Figures 7a and  b show the temporal evolution of the dislocation core dissociation angle, θ, and the core displacement in the e y -direction (P plane), ∆y, at T = 100 K under shear stress σ yz = 10 MPa.When the dislocation core has a P core (blue circles), ∆y increases quickly; i.e., P-glide is fast.On transformation to a π core (red circles), the dislocation pauses while "waiting" to transform back to a P core (Fig. 7a) upon which glide restarts on the P plane (Fig. 7b).Hence, the π core is "locked" (does not glide) and the P core is "unlocked" (glides easily).
The "locking" period predicted by kMC (∼ 10 −10 s) is much smaller than that observed in in situ TEM straining experiments (∼ 8 s) [10].There are two possible sources for this discrepancy.The TEM specimen is a thin foil, the free surfaces of which could provide strong drag on the dislocation.We have performed additional MD simulations, involving a dislocation line threaded at two free surfaces, to examine the drag effect (see SI for the simulation settings and results).The simulation results confirm that the free surfaces significantly reduce the dislocation mobility by imposing severe restrictions on core transitions.Second, our kMC simulations only study the motion of a short dislocation segment, rather than a long dislocation line.The differences in dislocation mobility between our work and TEM observations are attributed to the differences in dislocation line length.The motion of long dislocation lines involves collective motion of many dislocation segments, thus a higher free energy barrier should be overcome during the glide process.Moreover, long dislocations may migrate via kink pair nucleation and propagation, which necessitate inclusion of kink formation and migration energy barriers in estimation of dislocation mobility.Multiple core structures will likely be found on a long dislocation line.The interaction between these different core structures may also lower the dislocation mobility.This, coupled with the free surface restraint on core transitions (see SI) explains why short dislocation segments are more mobile than long dislocation lines in TEM observations.
To validate our kMC results, we simulated dislocation core motion at high temperature (500 K) at a high shear stress (σ yz = 60 MPa) by both kMC and DP-based MD (note that the MD time scale only allows the study of fast dynamics which can be achieved at high temperatures and high driving forces).Figures 7c and d show the evolution of θ and ∆y obtained through kMC simulations, while Fig. 7e and f show the same quantities under the same conditions from MD.The kMC and MD results are consistent.At this high temperature, the lockingunlocking mechanism is not easily seen, although it effectively lowers the core velocity along the P plane.

C. Dislocation Mobility
A shear stress σ = σ xz (e x ⊗ e z + e z ⊗ e x ) creates a PK force on the screw dislocation f = σ xz be x ; the dislocation moves, on average, in the e x -direction, i.e., on the B plane.We found that the dislocation velocity on the B plane, v B , is a linear function of σ xz (see the kMC and MD data in SI).The dislocation mobility on the B plane is M B = v B /(σ xz b).Alternatively, we may drive the dislocation motion on the P plane via shear stress σ = σ yz (e y ⊗ e z + e z ⊗ e y ) to obtain M P = v P /(σ yz b).If the energy barrier for a core transition or glide event is large, the intrinsic free energy barrier, Q i(κ)j , is the Peierls barrier and the dislocation motion is thermally activated.On the other hand, if it is small, Q i(κ)j cannot be interpreted as Peierls barrier; it is simply a parameter in the model which reproduces the frequencies obtained from MD. Dislocation mobility in the case of small energy barrier is phonon damping-controlled such that the viscous drag coefficient is M −1 B or M −1 P [42].Phonon damping is not explicitly modeled in kMC; rather it is captured through the parameterization of the frequency obtained from MD (similar to the method reported in Ref. [28]).
The temperature-dependencies of M B and M P are shown in Fig. 8. M B is much lower than M P at all temperatures, as suggested by Fig. 5f (where the dislocation core diffusion coefficient is a minimum/maximum at ϕ = 0/90 • .As shown in Fig. 8a and its inset, M B increases with temperature monotonically while M P exhibits a maximum at T = 300 K (indicated by the dashdotted line).The decrease of M P and increase of M B above 300 K is consistent with the MD results, i.e., the crosses in Fig. 8a inset.MD simulations from the literature [43][44][45][46] always show that the dislocation mobility decreases (or equivalently, the viscous drag coefficient increases) with increasing temperature.However, we note that MD results at low temperatures do not exist, since MD timescales do not suffice.The decrease in mobility (increase in viscous drag coefficient) at high temperature is usually interpreted as a phonon drag/damping effect.However, our kMC results suggest that glide on the P plane is also effectively damped by dislocation core transitions to B core.With increasing temperature, the B core is increasingly stable (Fig. 4a) such that the transition rate from the π core to the B core increases (Fig. 6b), leading to B-glide which contributes to zero motion on the P plane.We also investigated the effect of shear stress orientation and dislocation mobility anisotropy.A PK force may be applied in different directions, by choice of the relative magnitudes of σ xz and σ yz ; i.e., f = σ xz be x +σ yz be y .The orientation angle (maximum resolved shear stress plane, MRSSP) is χ = arctan(f P /f B ) = arctan(σ yz /σ xz ), where f P and f B are the PK forces resolved on to the P and B planes.We calculated M P and M B for various values of χ under 100 K and 200 K; the results are shown in Figs.7g and h.There is no surprise that M B decreases and M P increases as χ increases from 0 to 90 • , and the high-temperature mobilities are higher than the low-temperature counterparts.The dislocation mobility can be generalized to a tensor, M, defined by the relationship: v = Mf .The dislocation glide velocity component parallel to the PK force is where f = (cos χ, sin χ) T is the direction of PK force and the scalar dislocation mobility is where We measured the velocity component parallel to the PK force (v) and calculated the scalar mobility (M = v/f ) as a function of χ; see Fig. 8d.Next, we extracted M ij by fitting M vs. χ to Eq. ( 17); see the solid lines in Fig. 8d.We find that M 11 = M B (χ = 0), M 22 = M P (χ = 90 The glide probabilities are deduced from dislocation random walk (i.e., with no driving force), from which we extracted the intrinsic glide barriers.Applied stress does not change the intrinsic glide barriers, but rather biases the glide barriers according to Eq. ( 16).The Schmid effect is naturally included in this model, since the glide barrier will be lowered the most along the plane where the resolved shear stress is the maximum.
We do not investigate non-Schmid effects in this study as this is not an intrinsic property.To capture non-Schmid effects generally necessitates determination of how the entire 6-dimensional stress tensor alters the glide barrier.More specifically, it is possible to repeat our sampling and analysis presented in this work as a function of stress normal to the P-plane, π-plane or B-plane.This is beyond the scope of this study.

Long dislocation lines
The present simulations focused mainly on understanding the effect of intrinsic dislocation core properties/behavior on dislocation motion.Hence, all simulations were performed to short dislocation segments.While this helps obtain fundamental, intrinsic core behavior, it does not account for all aspects of dislocation dynamics.Nevertheless, it is of practical interest to understand how a long, screw dislocation line moves in α-Ti.The mechanism of the motion of a long dislocation line is associated with nucleation and propagation of kinks.Both kink nucleation and propagation necessarily involve the advance of local short dislocation segments.In this sense, the core properties extracted in this paper serve as essential input to such higher-level, long dislocation line models/simulations.
A serious treatment of long dislocations should be multiscale.Even the extant long (32b) dislocation MD simulations (e.g., see Ref. [14]) have not resolved the size effect issue.At a high temperature (k B T exceeds the kink energy), a dislocation line will likely undergo thermal roughening (i.e., the fluctuation amplitude of a dislocation line scales with the size of system) [47].If so, the size effect cannot be overcome by any finite length scale MD simulation.A possible strategy is to incorporate the intrinsic dislocation core properties as inputs for a multiscale method, such as kMC, rather than to simulate a long dislocation directly by MD.To do this, additional information is required (e.g., double kink formation and migration barriers or migration velocities); these may be obtained either directly from MD simulations or by derivation.For example, Edagawa's line tension model [48] provides a reasonable description of variations of the double-kink nucleation energy, which can also be directly obtained from atomistic simulation by modeling a kink structure [49].The kink migration barrier may be obtained from atomistic simulation for determination of Peierls barrier [49].If the kink migration barrier is much lower than the screw dislocation migration barrier, the kink velocity can be used in place of the migration barrier in kMC simulations [27].All above allow for the effective parameterization of long dislocations for kMC simulations (e.g., those proposed by Cai and collaborators [23,27]) of long dislocations in complex materials.

V. CONCLUSION
We have studied the finite-temperature core structures of ⟨a⟩ screw dislocations in HCP Ti, through a multi-scale framework.First, we characterize atomic interactions in Ti based upon machine learning, Deep Potentials (DP), which reproduce the stable/metastable dislocation core structures found via quantum mechanical, DFT calculations.DP was employed in molecular dynamics (MD) simulations of screw dislocation core structure at finite temperatures.MD provides the statistics of directional dislocation core dissociation (π, P and B cores) and directional core glide (π-, P-and B-glide).We found that the π core is stable and the P core is metastable, consistent with 0 K DFT results, while the B core is metastable above 300 K. Contrary to common understanding, the glide direction need not align with the core dissociation direction; e.g., π core can glide on the P plane.
The MD observations allow us to identify all important unit kinetic events associated with dislocation core motion.The events were categorized as either core transition (change in the dissociation direction) or core glide events (unit displacement along a slip plane).These events were incorporated into a kinetic model and that was parameterized through the MD data.The machine learning-based fitting procedure ensures that the frequency of each core structure and the translational and rotational diffusion coefficients produced by a kinetic Monte Carlo (kMC) simulation implementation of the model are consistent with MD data.We found that P core glide on the P plane event has the lowest core glide barrier, the transition between P and π cores has the lowest barrier among all core transition events, and the glide of the π core (on any plane) is very difficult.
With the parameters (barriers and frequencies) obtained by fitting, the proposed kMC simulation procedure is applicable to dislocation core dynamics at any temperature and applied stress.The dislocation will undergo a random walk (diffusion) in the absence of an applied stress; long-time dislocation core trajectories provide anisotropic translational and rotational diffusion coefficients.The former indicates that dislocation motion on the P plane is fastest and motion on the B plane is slowest.This implies that the π core is difficult to rotate and is stable while rotation away from the B core is fast.Under an applied stress, dislocation motion occurs through a locking-unlocking process at low temperatures, consistent with experimental observations.The locking behavior originates from the high energy barrier associated with π core glide.Application of different stress states yields that the motion of ⟨a⟩ screw core in Ti is anisotropic.The temperature dependence of this anisotropy is consistent with the (limited) MD predictions.
This work demonstrates that the intrinsic dynamic behavior of dislocations cannot be described based upon studies of 0 K dislocation core structures alone.Rather, statistical examination of the finite-temperature core structures is essential to determine, not only, the finitetemperature stability of different core structures, but also the kinetic properties of dislocation motion and core structure transitions.The kinetic model and parameters obtained in this work provide the necessary inputs for the higher-level approaches, such as kMC simulation of long dislocation line (for which kink nucleation and propagation are considered) and discrete dislocation dynamics.While the present study focuses on screw dislocation motion in Ti, the method described here is applicable to all dislocation types, at any temperature and stress state, in both simple (cubic) and complex (non-cubic) crystalline materials.

MD Simulations
The simulations employ a Deep Potential (DP) trained using DFT results for perfect crystals and defects in the α, β and ω phases of Ti [33].This DP successfully reproduces the 0 K core structures of the ⟨a⟩ screw dislocations in Ti as predicted by DFT [33].
The MD simulation cell geometry is shown in Fig. 1a.An HCP α-Ti single crystal cylinder is constructed such that the [1 210] is parallel to the cylinder axis; the Cartesian coordinate system employed has e x ∥ [10 10], e y ∥ [0001] and e z ∥ [1 210].The simulation cell is periodic along the e z -axis.An ⟨a⟩ screw dislocation, with both Burgers vector and line direction parallel to e z , is introduced in the center of the cylinder by displacing the atoms according to the anisotropic elasticity solution (lattice parameters and elastic constants for this potential as a function of temperature).The configuration is equilibrated at different temperatures in an N V T ensemble.The interactions between atoms in Region I (see Fig. 1a) are described by the DP for Ti.The atoms in Region II are described as an Einstein crystal, i.e., the atoms are tethered at the coordinates of the as-constructed (with anisotropic elastic displacements of atoms from their perfect crystal locations) configuration by harmonic springs, to avoid the free surface of Region II which will apply image force on the dislocation core.The interface between Region I and II has negligible effect on the simulation results.Details can be found in the SI.The spring constant is determined as 3k B T /⟨(∆r atom ) 2 ⟩, where k B is the Boltzmann constant, T is the absolute temperature and ⟨(∆r atom ) 2 ⟩ is the mean squared displacement of atoms at the temperature T .The radius of Region I is ∼ 160 Å, the width of Region II is ∼ 18 Å, and the dislocation line length is ∼ 6 Å.The whole system contains ∼ 33, 000 atoms; i.e., large enough that Region II and the interface between Region I and II have little influence on the random walk of the dislocation core about the center of Region I; see SI for details.The cylindrical sample containing a dislocation was equilibrated for 2 ns at each temperature.All MD simulations were performed using lammps [50].
The dislocation configurations were thermally equilibrated at temperatures 300-900 K (well below the HCP-BCC transition temperature, 1250 K).The atomic configuration was recorded every 50 fs.Energy minimization for 10-20 steps (by conjugate gradient with line-search step size 0.01 Å) was conducted for each of these atomic configurations to remove thermal vibration in order to clearly visualize/analyze the atomic structure [23].Such energy minimization has little effect on the dislocation core distribution (for details, see SI).

Nye tensor parameters
The Nye tensors are visualized with perfect crystals as references.Two key parameters govern the presentation of the Nye tensor plot are the cutoff distance for constructing a neighbor list and the maximum angle (Θ) employed to identify matches between p and q vectorshere, p and q denote the radial distance vectors between each atom and its neighbors in the reference and current systems, respectively.In our investigation, we have chosen a cutoff distance that corresponds to 1.3 times the equilibrium lattice constant at the temperature of interest.We set Θ to 10 • .

Prediction of event frequencies associated with dislocation dynamics
In S1 of Sect.III C, we need an initial guess of the frequencies of all events at all temperatures based on the MD results.We treat the core transition events (i ̸ = j) and the core glide events (i = j) differently.
The frequency for a core transition event i(κ)j (i ̸ = j), ν i(κ)j , is obtained from the MD results by where ∆t is the time interval between recording atomic configurations, N i is the number of configurations showing the i core, N i(κ)j is the number of the i(κ)j events recorded, and P i(κ)j is the conditional probability P(i → j|i) = P(i → j)/P(i) = N i(κ)j /N i .The sampling interval ∆t is chosen small enough that few events are missed but not so small that thermal vibrations lead to incorrect event identification.We set ∆t = 0.1 ps (i.e., close to the inverse Debye frequency for Ti).
The frequency for a core glide event i(κ)i, ν i(κ)i , may also be extracted from the MD simulation atomic configurations, sampled with time interval ∆t.Suppose that the number of the i(κ)i events is N i(κ)i and the displacement vector of the m th i(κ)i event is d i(κ)i (m).Then, the total time spent on the i(κ)i-glide is N i(κ)i ∆t and the glide distance of the i(κ)i event is s κ • d i(κ)i (m), where s κ = n κ × b/|b| (n κ is the normal to the κ plane and b is the screw dislocation Burgers vector).The frequency for the i(κ)i-glide event is thus obtained from where L i(κ)i is the shortest glide distance, i.e., one period of Peierls barrier along the direction s κ .s B , s P and s π are the unit vectors in the [10 10], [0001] and [10 12] directions, respectively.According to Figs. 3a-f, L B(B)B = √ 3a/2, L P(P)P = c/2, L π(P)π = c and L π(π)π = 3a 2 /4 + c 2 .

3 e
FIG. 1. MD simulations of the ⟨a⟩ dislocation core structure in α-Ti.(a) The cylindrical model and Cartesian coordinate system employed in the MD simulations.The blue "+" symbol indicates where the screw dislocation was introduced.The atomic interaction in Region I and Region II are all modeled using DP, while atoms in Region II are treated via an Einstein model.(b) An HCP unit cell, where the P, π, π ′ and B planes are identified.(c)-(g) The dislocation cores observed in the MD simulations (following 10-20 energy minimization steps -for improved visualization, see Method and Ref. [23]).The dark and light gray circles denote atoms on successive (1 210) planes (perpendicular to the ez-axis) in the perfect crystal.The arrows represent the differential displacement map while the color heat map shows the αzz distribution.Blue "+" symbols indicate the core positions determined by the first moment of αzz.
FIG.2.Dislocation core trajectory, dissociation direction and glide direction.(a) The example trajectory of a dislocation core obtained by the MD simulation at 500 K for 2 ns.The blue, red and green points denote the positions where the core structures are P, π and B, respectively.The inset shows that the P/B cores and the π cores should be distributed on alternating P planes.(b) Temporal evolution of the core dissociation direction θ extracted from the same MD simulation.(c) Distribution of θ at the temperature ranging from 300 K to 900 K.The horizontal lines denote the boundaries we defined to distinguish the core structures at 500 K. (d)-(f) Distributions of the deviation of the glide direction ϕ away from the dissociation direction θ for the B, π and P cores.(g)-(i) The differences of the αzz maps at two moments, corresponding to the four dislocation core glide events observed in the simulation.The blue/red "+" symbol indicates the core position in the previous/next moment.

FIG. 4 .
FIG.4.The equilibrium probabilities of core structures and the free energy differences.(a) The probabilities of the π (red circles), P (blue squares) and B (green triangles) cores obtained by the MD simulation (open symbols) and the kMC simulation (solid symbols).(b) The core energy differences, ∆FPπ (blue squares) and ∆FBπ (green triangles), obtained from the MD data in (a).The solid lines are the fitting results based on the formula Eq. (5).

FIG. 5 .
FIG. 5. Mean squared displacement and diffusion coefficients.(a) and (b) show the mean squared displacements (MSDs) of a dislocation core in the x-and y-directions at different temperatures.(c), (d) and (e) show the mean squared angular displacements (MSADs) of a core dissociation direction starting from θP, θπ and θB.(f) Translational diffusion coefficient as a function of glide direction, D T (ϕ) at different temperatures.(g) Rotational diffusion coefficient as a function of dissociation direction, D R (θ) at different temperatures.In all figures, the circles denote the MD data, the solid lines denote the results of kMC simulations based on the optimized frequencies {ν i(κ)j }, and the dotted lines denote the results of kMC simulations based on the fitted parameters {ν 0 i(κ)j } and {Q 0 i(κ)j }.

− 3 K − 1 )FIG. 6 .
FIG.6.Kinetic model and coefficients.(a) Schematic of the energy landscape for the i(κ)j event.Fi and Fj are the free energies of the i and j cores; ∆Fij = Fj − Fi and ∆Fij = 0. F b i(κ)j is the total free energy barrier and Q i(κ)j is the intrinsic barrier for i(κ)j.(b) η i(κ)j (Eq.(12)) vs. inverse temperature.(c) η i(κ)i (Eq.(12) for glide events) vs. the inverse temperature.In (b) and (c), symbols and solid lines denote MD data and fitting results (Eqs.(12) and (13)).

FIG. 7 .
FIG. 7. Temporal evolutions of the dislocation core dissociation angle (θ) and the core displacement along P plane (∆y).(a) and (b) The core trajectories under σyz =10 MPa at 100 K predicted by kMC.(c) and (d) The core trajectories under σyz =60 MPa at 500 K predicted by kMC.(e) and (f) The core trajectories under σyz =60 MPa at 500 K obtained by MD.In all figures, blue, red and green circles label the P, π and B cores, respectively.

8 FIG. 8 .
FIG. 8. (a) Arrhenius plot of dislocation mobilities on the P plane for σyz (blue) and B plane for σxz (red).The circles and crosses are, respectively, the kMC and MD data.The dashdotted line indicates T = 300 K. (c) and (d) The mobilities on the B and P planes (MB and MP) as a function of PK force direction χ.(d) The scalar dislocation mobility as a function of χ.The symbols are kMC data while the solid lines are fits to Eq. (17).In (b)-(d) the blue and red circles correspond to 100 K and 200 K.
The second moment of α zz is the tensor [C ij ], where C 11 = x 2 , C 22 = y 2 and C 12 = C 21 = xy.Suppose that the (normalized) eigenvector of [C ij ] corresponding to the maximum eigenvalue is e * .Thus, the core dissociation direction at t is given by θ(t) = arctan e * y (t) e * x (t) .
[38]nd M 12 is negligibly small in comparison with M 11 and M 22 .Since M 11 ̸ = M 22 , the dislocation velocity v is, in general, not in the same direction as the PK force f .The off-diagonal component M 12 relates to how the PK force, resolved on B/P plane, influences dislocation glide on the P/B plane -this is a non-Schmid effect[38].M 12 ≈ 0 indicates that dislocation glide is well-described by the Schmid law in our kMC model.