Analyzing the cross slip motion of screw dislocations at finite temperatures in body-centered-cubic metals: molecular statics and dynamics studies

The plasticity of body-centered-cubic metals at low temperatures is substantially determined by the screw-dislocation kinetics. Because the core of screw dislocations in these metals has a non-planar structure, its motion is complex. For example, although density functional theory predicts slip on a {110} plane, the actual slip plane at elevated temperatures differs from the prediction. In this work, we explored state-of-the-art atomistic modeling methods and successfully reproduced the transition of the slip plane through a temperature increase. We then devised an algorithm to analyze the activation of dislocation jump over the Peierls barrier and discovered a possible origin of this unexpected phenomenon: thermal fluctuation leads to the kink-pair nucleation for cross slip jumps with no transition of the dislocation core structure.


Introduction
Transition metals and alloys, such as steels, are ubiquitous in the infrastructure of automobiles, ships, bridges and nuclear power stations. Accordingly, their mechanical degradation may have catastrophic consequences. It is thus essential to understand the plasticity of these materials, which is mainly governed by the dislocation motion. Therefore, having a complete picture of the physics of the dislocation motion is fundamental to understand their mechanical degradation.
In body-centered-cubic (BCC) metals, the plastic deformation at low temperatures is accompanied by screw-dislocation glides whose motion occasionally seems complex [1][2][3]. The non-planer structure of the dislocation core is believed to be responsible for such complexity [4][5][6]. The threshold barrier for a screw dislocation jump, also known as the Peierls barrier, is so high in BCC metals that the dislocation line is thermally activated from a Peierls energy valley to shift to the next via the kink-pair (or double-kink) nucleation and kink migration [7][8][9]. This kinetics is understood as a stochastic process and is supposedly represented by an Arrhenius-type rate equation following the transition state theory [10]: , where γ is the dislocation velocity; γ 0 is a constant; ΔH the Peierls barrier; k the Boltzmann factor; and T the absolute temperature. This stochastic theory is supported by the fact that the brittle-to-ductile transition, regulated by the activation of the screw dislocation motion, follows a stochastic theory [11].
Even though this theory has been proposed, its detailed understanding remains physically unclear. For example, some experiments reveal that the screw dislocations at low temperatures have a jerky motion [12] that cannot be explained by the stochastic theory. Moreover, the slip plane switches from {110} to {112} as the temperature increases [13]. For example, Francioci et al [14] experimentally found that the slip plane in pure iron (Fe) can be straight on a {112} plane. In addition, Marichal et al observed 'sharp slip' traces on the {112} plane [15]. Seeger et al [16] claimed that they succeeded in explaining the presence of the pure {112} plane. However, afterwards the {112} glide plane proved to be an average plane composed of two different {110} glide planes [15,[17][18][19], so that the glide occurs in zig-zag order: 101 110 101 110 101 110      ( ) ( ) ( ) ( ) ( ) ( ), refuting the theory of Seeger et al [16]. Until now, the cause of the appearance of the {112} plane has not been explained.
Even with the state-of-the-art technology, there are no experimental means to further analyze and obtain the atomistic understanding of experimental observations, such as the appearance of the {112} plane, because the dislocation trace observations have neither atomic-scale spatial precision nor picoseconds time resolution to observe atomistic events.
To explore the fundamental complexity mentioned above, atomistic-scale computational methods are required. First principles calculations based on the density functional theory (DFT) are applied to studies on the stability and mobility of screw dislocations [20]. For example, some recent reports reveal not only the one-dimensional Peierls barrier of tantalum [21] and other BCC metals [22], but also the two-dimensional (2D) Peierls potential in BCC iron [18], tungsten [23], and other BCC transition metals [24]. Different studies predict the kink-pair formation energy [24,25] and the trajectory [26] of screw dislocations in BCC transition metals. Despite these achievements, the DFT method is practically restricted to molecular statics and not directly applicable to tracking mobile defects at finite temperatures, which are inevitably found when analyzing the plasticity of materials. Attempts to include magnetism into empirical potentials deserve attention [27][28][29][30], because they succeeded in predicting correctly the core structure and the energy barrier of screw dislocation migration. Nevertheless, they are not applicable to tracking a decent length of dislocation at finite temperatures. Accordingly, it is necessary to resort to molecular dynamics simulations utilizing simpler empirical inter-atomic potentials, e.g. based on the embedded-atom method (EAM). So far, this classical method has been successfully applied to the observation of a screw-dislocation motion activated by the kink-pair nucleation [31,32]. Furthermore, the nudged-elastic-band method allows for a more detailed observation [33]. In addition, the dislocation mobility [17] and kink migration [9] have also been discussed using empirical potentials. Regarding the appearance of slips on {112}, Chaussidon et al [34] claim that the degenerated easy-core configuration causes this slip plane but some studies deny this notion, indicating that the {112} slip plane is related to the metastability of the split-core configuration [35,36]. Thus far, studies based on the empirical potential method analyze the dislocation motions with a spatial resolution as large as the Peierls valley interval. Accordingly, it cannot be observed how a kink pair at finite temperatures is dynamically nucleated, and this inaccuracy restricts the understanding of the cause of cross slips necessary for the {112} slip plane observed at elevated temperatures. To overcome this issue, it is necessary to be able to visualize critical moments with a much higher spatiotemporal resolution than the previous studies.
In this study, we successfully devised a numerical method to track the 1/2〈111〉 screw dislocation kinetics in α-Fe with unprecedentedly high spatiotemporal resolution, and analyzed the nucleation of kink pairs under a constant shear strain. The results indicate that the dislocation glide motion on {110} is dominantly observed at low temperatures, agreeing with experimental observations. The results also show the {112} slip plane at high temperatures, also in accordance with experimental facts. Our dislocation tracking algorithm demonstrates that the thermal fluctuations of the dislocation line causes kink-pair nucleation for cross slip jumps and triggers the departure from the {110} slip plane. This suggests that the {112} slip plane is not due to a change in the core structure beyond a critical temperature, but rather to thermal fluctuations of the dislocation line. The results of this computational study also show an intriguing property of dislocation motions: a single screw-dislocation core structure can cause two different slip planes depending on the temperature.

Computational methodology
We conducted two types of atomistic simulations of α-Fe crystals using the molecular dynamics simulation package LAMMPS [37], i.e. molecular statics to analyze the 2D Peierls potential for 1/2〈111〉 screw dislocations, and molecular dynamics for the dislocation motions. In both cases, a BCC Fe crystal was created in the simulation box (figure 1), where the 〈111〉 direction was parallel to the x-axis and an anti-clockwise 1/2〈111〉 screw dislocation was inserted along this direction. Free surfaces were adopted in the z-direction, and the remaining boundary conditions were decided case by case. For the simulation box, typical x-, y-and z-side lengths were set to 195, 322, and 159 Å, respectively. Occasionally, we changed the box dimension to verify the interpretation of certain numerical results.
At finite temperatures, the dislocation core trajectory on the {111} plane may deviate from its ideal track due to thermal disturbances, and the 2D Peierls potential of 1/2〈111〉 screw dislocations becomes critical to analyze this thermalized trajectories. The 2D Peierls potential for the adopted empirical inter-atomic potential was calculated as follows: we inserted an approximate screw dislocation normal to a {111} plane, i.e. normal to the x-plane in the simulation box (figure 1), into the perfect crystal by artificially displacing some atomic blocks in the x-direction; then, the whole simulation box was relaxed to allow for a physically stable linear dislocation at an easy-core configuration and the total energy was then measured. Next, we artificially shifted the atoms near the core in the x-direction and changed the core position to anther configuration; the system was then relaxed with fixing atoms for the dislocation core vicinity, and the total energy was again measured. The same procedures were iterated to obtain enough energy values to draw the global 2D potential.
For the screw-dislocation motion analyses, we performed deformation molecular dynamics simulations with a constant shear strain after inserting a linear screw-dislocation and equilibrating the system at a given temperature. As depicted in figure 1, we created two sub-blocks at the top and bottom of the simulation box: the bottom block is always fixed at the original position, while the top one is displaced with a constant rate toward the positive xdirection, i.e. [111], and the block between the top and bottom is gradually deformed. For all the cases, we conducted an integration with constant number, volume, and energy (namely, a constant NVE integration) to update the atomistic configuration with a 1 fs time step interval.
High-accuracy tracking of screw-dislocation core positions is the central issue of this study. For this purpose, we devised the following algorithm. At the beginning of each simulation, before inserting a screw dislocation, we decompose the perfect crystal into many atomic blocks, each of which is a thin board parallel to the x-plane with the x-thickness of b, the length of Burgers vector, so that the atoms in each board cover the whole x-plane of the simulation box. Then, the positions of all the atoms in each thin board are recorded. After inserting a linear screw dislocation into the perfect crystal by using the method mentioned above, we measure the x-displacements u x for all the atoms in each thin board based on the recorded original positions. We collect u x i , where i is the index of an atom. We can then obtain the differential displacement map of a thin board. In this map, we search for three triangle lattice points with non-zero dislocation helicity; this triangle approximately locates the dislocation core position of the thin board. The exact position of a core center inside the triangle can be calculated from the differential displacements of the three lattice points using equation (2) in [18], if the system is fully relaxed. We, however must track the mobile dislocation motion, and this method sometimes fails when a core center is located at the border of two neighboring triangles. To avoid such discontinuity, we slightly revised the method, i.e. we determined the core center position using the differential displacement map of a wider area. The details of this revised method are reported elsewhere [18]. By connecting the core positions for different thin boards, we obtained the whole view of a dislocation line. The molecular dynamics simulations create not only linear dislocations but also curved ones. Obviously, the above method can be applied to the curved ones because each dislocation segment position can be determined independently. In addition, we can trace a dislocation motion by accumulating consecutive snapshots of the dislocation line. This method works at low temperatures, but it sometimes yields erroneous results as the temperature increases because lattice vibrations blur the dislocation core structure. To avoid this, we measured u i (t) using the time averages of the atomic positions instead of a single snapshot. This movingaverage technique was adopted for all the cases, and the the 10 and 1000 fs time window was determined by try-and-error.

2D Peierls potential analyses
We chose pure crystalline α-Fe as the BCC metallic system, because, among all the BCC metals, α-Fe has been the most thoroughly investigated using molecular dynamics, e.g. [17,[31][32][33]. Choosing a good empirical potential is essential to atomistic simulations because many α-Fe potentials fail to represent the correct 1/2〈111〉 screw dislocation. The EAM potential called Mendelev Type-2 potential [38] has been used in some previous studies [17,31,32,34] as it correctly reproduces the compact core (or the non-degenerated core) structure when the system is fully relaxed. We consider that the most important criterion is that the potential must correctly represent the single-hump Peierls barriers for a single 1/2〈111〉 screw dislocation jump, as predicted by the DFT results [18]. Otherwise, the screw dislocation in this system is expected to rest at an intermediate configuration during a single jump, which is physically incorrect. The high resolution analysis of screw-dislocation motions is the main interest of this investigation, and wrong Peierls barriers can easily damage the whole study. According to figure 5 of Ventelon et al [23], which compares Peierls barriers of the screw dislocation for different EAM potentials, only a potential called MCM2011 [25] has a single hump, for which we anticipate the same kinetics as predicted by DFT. To our best knowledge, no molecular dynamics study on the screw dislocation motion using this potential has been reported. To further analyze the EAM potential, we visualized its 2D Peierls potential for a 1/2〈111〉 screw dislocation on the {111} plane normal to the dislocation line using the method described in the previous section. Such analyses are not new as the same results are given in [24], but we revisit this problem to clarify merits and demerits caused by adopting the EAM potential. Figure 2(a) shows the 2D Peierls potential landscape predicted by the DFT results under no external stress. To obtain this diagram, we transformed the DFT evaluations of the screw dislocation reported by Itakura et al [18] inside the minimum triangle area into a larger area. The global minima of this landscape are located at the easy-core configuration ('E' in figure 2(a)); if no stress is applied at extremely low temperatures, the screw dislocation core is supposed to be stabilized at the center of the triangles composed of three neighboring lattice points on the {111} plane. There is another type of triangles whose center points are called the hard-core configuration ('H' in figure 2(a)), at which the corresponding Peierls potential has an inconspicuous local maximum. The split core ('S' in figure 2(a)) practically located almost at the lattice points on the {111} plane has the maximum Peierls energy. Thus, the screw-dislocation core is likely to travel from an easycore configuration to the next on the {111} plane as schematically shown in figure 2(a). Additionally, figure 2(b) shows the 2D landscape evaluated for the MCM2011 EAM potential. Similar to the DFT prediction, its global minima are located at the easy-core configuration ('E' in figure 2(b)). Here, the main difference from the DFT prediction is that the hard-core configuration has the maximum energy, while the split-core configuration corresponds to an inconspicuous local maximum. Thus, we must expect that the screw dislocation migration is qualitatively different from what we expect from the 2D Peierls potential given by the DFT . The optimum migration path is expected to pass above the splitcore configuration, as represented by the dashed line in figure 2(b). Moreover, the path is significantly curved, while for the DFT prediction it is approximately linear. A critical question is whether or not these differences mentioned above render the EAM potential useless for our purpose. Overall, the split-and hard-core configurations in the EAM potential function as the hard-and split-core configurations in the DFT prediction, and the dislocation core travel on a 2D Peierls potential is similar to that of the DFT prediction. Regarding the curved path for the MCM2011 potential, Dezerald [26] reported that the DFT-predicted path of niobium (Nb) is also curved, and the slip-plane transition for Nb crystals occurs above 175 K, similar to other BCC metals [13]. This means that a curved path does not necessarily affect our purpose. Accordingly, we expect that the screw dislocation modeled by the MCM2011 EAM potential travels on a 2D Peierls potential landscape that satisfies the inherent properties of BCC transition metals. For comparison, 2D Peierls potential for the Mendelev Type-2 potential is shown in figure 2(c). Its general landscape is similar to that for MCM2011, i.e. the minimum and maximum dislocation energy values are respectively at the easy-and hard-core configurations. Different from the MCM2011 potential, the split-core configuration has the local minimum and the optimal migration trajectory is expected to touch the split-core configuration as illustrated by the dashed line in figure 2(c). This local minimum also causes a double-hump Peierls barrier [23], and the dislocation line is likely to rest at this intermediate configuration during the single jump, which is essentially different from what we expect from the DFT results. Additionally, the stable split core promotes the dislocation to attract the trajectory to the lower layer of the glide plane, i.e. it promotes the cross slips and causes the slip-plane transition. Therefore, the Mendelev Type-2 EAM potential has the inherent property of promoting {112} slip plane, and may not be suitable for examining the {112} slip plane's origin, which is our prime interest. In fact, even if the maximum resolved shear stress plane is set at {110}, molecular dynamics simulations with this potential predict {112} slip planes [31,34], with being a single exceptional study [32], in which the {110} slip plane is predicted.

Molecular dynamics results for dislocation mobility
The parameter domain of the selected molecular dynamics simulations are shown in figure 3: we chose four temperatures, 20, 100, 200 and 300 K; and five strain rates from 1.67×10 6 to 1.67×10 10 s −1 for our principal study. Among these only limited case studies were performed for 200 K as seen in figure 3. 20 K is extremely low, but it is convenient to observe dislocation motion because thermal noise hardly disturbs the observation. For the typical cases, we set the periodic and free boundary conditions along the x-and y-directions, respectively. As seen in figure 3, the basic dislocation behaviors for the fastest deformation (1.67×10 10 s −1 ) were different from the four slower cases (1.67×10 6 to 1.67×10 9 s −1 ) at the same temperature. Thus, we believe that the fastest deformation causes some artifact. Contrastingly, we did not observe a significant deformation-rate dependence on the basic dislocation behaviors among the slower rates, hence they were further analyzed. As indicated in figure 3,  We used the post-processing algorithm explained above to examine the dislocation migration. Figures 4(a) and (b) show the dislocation behaviors at low temperatures, in which the slip plane is {110}. At 300 K, the slip plane departs from {110} and the dislocation line becomes jagged because of the increased thermal perturbation (see figure 4(c)). It is worth noting that the cross slips occasionally occur even at 20 K as rare events. The core trajectories on the {111} plane for 100 and 300 K are compared with each other in figure 5, where the slip-plane transition from {110} to {112} as the temperature rises is clearly observed. Despite the differences between the adopted empirical potential and the DFT prediction (figure 2), the molecular dynamics simulation succeeded in reproducing the experimentally observed slipplane transition [14,15]. Based on additional sensitivity analyses, we confirmed that this transition occurred irrespective of the adopted dislocation lengths, i.e. the length of the x-side of the simulation box. Figure 5 also indicates that the {112} glide is a wavy trace composed of two different {110} slips that randomly appear; agreeing with previous theoretical and experimental studies [15,[17][18][19].
As previously mentioned, the prime motivation of this study is to find the origins of the {112} slip plane at elevated temperatures. Therefore, we analyzed the moment of the cross slip. Figure 6 shows the dislocation development a high time resolution, suggesting that the fluctuation of the dislocation line triggers a kink-pair nucleation for a cross slip: part of the dislocation largely fluctuated backward from the remaining part at t=0 ps (figure 6(a)) triggering the cross slip kink pair. As already mentioned, the cross slips occasionally occur at  4(a)); this clearly depicts the moment of a cross slip due to low thermal noise ( figure 7). In this diagram, a cross slip motion also starts from a backward fluctuation (see the diagram at t = 0 ps in figure 7(a)); then causing the nucleation of the cross slip kink pair at t=0.6 ps ( figure 7(b)), which is followed by kink migration at t=1.2-3.6 ps (figures 7(c), (d)). This behavior suggests that a large fluctuation of a dislocation part is the origin of the cross slip motion. The fluctuation is enhanced at higher temperatures, thus the frequency of   In addition to the above results that correspond to dislocation motions with the periodic boundary condition along the x-direction, it is important to evaluate the influence of the boundary on the slip-plane transition. For this, we performed additional simulation studies with free boundary conditions and confirmed that the slip-plane transition from {110} to {112} also occurred as the temperature increases from 100 to 300 K ( figure 8). Accordingly, we did not observe any significant differences between the periodic and non-periodic boundary conditions in the x-direction. Figure 9 shows snapshots of the dislocation when the cross slip occurs at 20 K, where the thermal noise is suppressed, with free boundary conditions along the x-direction. As seen in figure 9, we found that the kink pairs triggering the cross slips occurred when part of the dislocation was largely delayed from the rest. This further supports the fact that thermal fluctuation causes the cross slip motion thus resulting in a {112} slip plane.

Core trajectory on the {111} plane
In the previous section, we showed that the fluctuating part of the dislocation line causes a kink-pair nucleation toward the cross slip direction. This behavior can be analyzed from a different point of view, i.e. from the trajectory of the dislocation core on the {111} plane normal to the Burgers vector. The trajectory at the {111} plane, where a kink pair for the cross slip is nucleated in figure 6, is shown in figure 10, in which the two fundamental slips on 110 ( ) and 011 ( ) can be compared. Without the fluctuation, the dislocation core rests around an easy-core configuration, then gradually approaches the optimal path above the split-core denoted by blue circles, and the 110 ( ) slip occurs. However, the fluctuation causes the core position to pass through the easy-core configuration, and a kink pair for the cross slip is nucleated. From this, we confirm that the fluctuation leads to the kink-pair nucleation for the cross slip. We thus consider that the modest peak of the 2D Peierls potential at the split-core configuration promotes the cross slip. The cross slip is not related to the transition of the stable core structure as the temperature increases because it occurs even at 20 K.  It is also important to consider how the selection of empirical potentials influences the molecular dynamics simulation. For this, we conducted additional studies with the Mendelev Type-2 potential [38] that is adopted in most of the previous studies [17,31,32]. In general, the results were similar to those given by the MCM2011 potential, i.e. the slip plane is {110} at 20 K, and it shifts to {112} at 300 K. Also, the probability of the cross slip slightly increases at 100 K. This change is anticipated because of the 2D Peierls potential around the split-core configuration (see figure 2(c)), i.e. the 2D Peierls potential of the Mendelev Type-2 EAM potential has no peak at the split-core configuration. Therefore, it can be said that the the results of this study are independent of the fine tuning of the empirical potential as there are only slight quantitative differences. On the contrary, we expect less cross slip probabilities for the DFT 2D Peierls potential, i.e. a higher slip-plane transition temperature is expected because its 2D Peierls potential has a straight optimal path as depicted in figure 2(a).

Influence from the boundary conditions of the simulation block
The boundary conditions in the x-directions were already discussed; yet the other directions. For the y-direction, i.e. the dislocation progressing direction, previous molecular dynamics studies [17,31] adopt the periodic boundary condition. In this study, we adopted the free boundary condition, which we believe to be physically natural as the simulation box had a yside as large as 230 Å, and the dislocations away from the both boundaries were analyzed. However, it is interesting to evaluate how this periodic boundary condition influences the dislocation behavior. Several analyses show that the dislocation glided on {110} over the temperature range of 20-300K; and the cross slip was significantly suppressed by the periodic boundary condition along the y-direction. Accordingly, this condition rather restricts the freedom of the dislocation glide motion. We observed cross slips at 600 K, but they occurred both ways, maintaining the average slip plane on {110}; hence the periodic boundary condition did not reproduce the experimental observations.
To impose shear stress to the dislocations, we fixed the upper and lower boundaries in the z-direction as shown in figure 1, but this rigidity may create some artifacts. To reduce this influence, we doubled the z-side of the simulation box in additional studies. We found that the Figure 10. The dislocation core path on a {111} plane at 300 K and the strain rate of 6.68×10 9 s −1 . A path of a dislocation fragment that caused the kink-pair indicated in figure 6 is shown; the red triangles are plotted every 100 fs, and the black triangles are plotted every 1 ps; blue circles are the lattice points; the arrow shows the trajectory when the fluctuation triggers the kink pair for the cross slip.
slip-plane transition between {110} and {112} was still observed as the temperature increased; demonstrating that the slip-plane transition is irrelevant to this boundary condition.

Conclusions
The kinetics of the 1/2〈111〉 screw dislocation in α-Fe were computationally simulated under constant shear strain conditions using molecular dynamics. The results indicate that the dislocation dominantly glides on the {110} plane at low temperatures. Moreover, as evidenced in the many experiments, the slip-plane transition to {112} was observed as the temperature increases. We also devised a numerical method to follow the screw dislocation with unprecedentedly high spatiotemporal resolution, and analyzed the dislocation behavior given by the molecular dynamics simulation. Our method demonstrated that the thermal fluctuation of dislocation lines causes kink-pair nucleation for the cross slip and triggers the departure from {110} slip plane, suggesting that the slip-plane transition is not due to the transition of the core structure over a critical temperature. This means that the identical dislocation core with the same maximum resolved shear stress plane can cause different slip planes depending on the temperature.