Our MD simulations were performed using the graphics processing unit (GPU)-based MD package MDPSCU^[28] (with the sources available at http://radphys.scu.edu.cn:8088). We used the many-body semi-empirical potential of the Finnis-Sinclair type proposed by Ackland et al.^[26] for Zr–Zr potential. This Zr–Zr potential has been widely used to study the stable configurations and diffusion of SIAs in Zr.^[3,10,11,14] Although the mostly stable SIA configuration predicted by the potential is the basal-crowdion (BC), which deviates from the basal-octahedral (BO) predicted by the most recent ab initio modelling, the prediction that the basal-split (BS) configuration has the second-lowest energy agrees with the ab initio prediction.^[15–17,21] The potential used in the present paper makes predictions comprehensively closer to those of the newest ab initio modelling than a later-developed potential.^[24]

The initial configurations of the simulation boxes were prepared by introducing a Zr atom into a Zr substrate in an hcp structure at a randomly selected position. For statistical purposes and optimal use of GPU-accelerated MD simulations,^[28] we generated 1000 independent replicas of simulation boxes, each of which contained 6481 Zr atoms with a size of 15a₀ × 15.59a₀ × 19.14a₀, where a₀ represents the lattice constant of Zr. The x, y, and z axes were set along , , and , and the periodic condition was applied to the three directions. The boxes were then thermalized and relaxed at a given temperature ranging from 300 K to 1200 K until thermal equilibrium was reached. The boxes were relaxed further for 30 ps in time-steps of 0.5 fs. The instant configurations of the simulation boxes were generally recorded every 0.1 ps. A finer recording time step was used only for observing the detailed transition path between SIA configurations. Because 1000 replicas of simulation boxes were used in parallel, the simulation procedure is equivalent to running one simulation box for 30 ns. The time scale is probably not long enough for calculating the diffusion coefficient with using the Einstein formula but is sufficient for calculating the frequencies of jump events defined below.

Methodologically, correct identification of SIAs is essential for a quantitative analysis of SIA migration features, particularly in high temperature cases. One approach to identifying SIA is to search spherical regions of a given cut-off radius centred on perfect hcp lattice sites. An atom is identified as an interstitial atom if it is not within any of the spherical regions. The shortcoming of this approach is that the search results may be sensitive to the cut-off radius, especially at high temperatures. For example, multiple SIAs could be identified even if only one SIA were introduced into a simulation box. Thus, we used an alternative approach in which the Wigner–Seitz cells of the hcp lattice sites are searched. The Wigner–Seitz cells are the Voronoi polyhedrons centred on the lattice sites. By the definition of Voronoi polyhedron, an atom is within the Wigner–Seitz cell of a lattice site when the atom is closer to the lattice site than to any other lattice site. Using the linked-cell technique and GPU parallel computing, the search process can be quickly completed. Because the Wigner–Seitz cells are space filling, there must be at least one Wigner–Seitz cell containing more than one atom if an SIA is introduced into the simulation box. We defined the Wigner–Seitz cell containing more than one atom as the location of self-interstitial atom (LSIA). This approach was used to identify interstitials and vacancies produced in cascade collisions and was shown to be robust.^[29] We also tested the robustness of the approach to identifying LSIA. In the simulation box created above, only one LSIA could be found, and an LSIA contained only two atoms for all temperatures considered in the present paper. With respect to multiscale simulations, in which the LSIA extent is very smaller than the distance that an LSIA moves before its encounter with other defects, it is not important to distinguish which atom of the two atoms at the LSIA is the self-interstitial atom. Thus, instead of identifying the detailed atomic events occurring during SIA migration, we use the change of the LSIA to describe the SIA migration.

Although the LSIA migration is a coarse grained description of the SIA migration, the migration of the LSIA can be anisotropic due to the dumb-bell like structure along direction formed by the two atoms at the LSIA. To quantifying the anisotropic LSIA migration, we defined three types of LSIA jumps as schematically displayed in Fig. 1. As mentioned above, the mostly stable SIA configuration predicted by the potential of Ackland et al.^[26] is the BC configuration. It is observed in our MD simulations that for most of time, the LSIA performs frequent back-and-forth swinging between a pair of neighbouring lattice sites in the same basal plane, for example the L1 and L2 site (Fig. (1a)). For the convenience of description, we call such a pair of lattice sites the LSIA pair. It would overestimate the number of jump events if we simply define the change of LSIA as a jump. Thus, we define the change of LSIA pair as an LSIA jump. For example, the lattice site pair L1–L2 is an LSIA pair at first (Fig. (1a)). With the passage of time, there is a probability that the LSIA pair could change from L1–L2 to L2–L3 (Fig. (1b)), This change indicates the occurrence of an LSIA jump. Such a jump that remains in the same direction of the immediately previous LSIA pair (L1 to L2) is referred to as an in-line jump (ILJ). It is also possible that the LSIA pair could change from site L1–L2 to L2–L4 (Fig. (1c)). In this case, the direction of LSIA pair L2–L4, [ , differs from the immediately previous direction (from L1 to L2). Such a jump involving a directional change from one direction to another equivalent direction on the same basal plane is defined as an off-line jump (OLJ). The change of LSIA from site L2 to L5 may also possibly occur (Fig. (1d)), which is a jump from one basal plane to another defined as an off-plane jump (OPJ). The robustness of calculating the number of jump events was also checked by that the same number of jump events is obtained when different time intervals, 0.1ps and 0.01ps, are used to detect the jump events.

Fig. 1.

Figure Option ViewDownloadNew Window

Fig. 1. (color online) Schematic graphs for the movements of interstitial atoms and LSIA. Grey spheres denote the atoms at normal hcp lattice sites at which their Wigner–Seitz cells contain one atom. Dots, for example site L1 in panel (a) and site L2 in panel (b), denote LSIAs at which two atoms (red spheres) fall in their Wigner–Seitz cells. The normal lattice sites nearest to an LSIA are LSIA candidates.

According to the simulation approach described in Subsection 2.1, MD simulations were performed at different values of substrate temperature T_s. As an example from the 1000 simulation boxes, figure 2 displays the typical trajectories for T_s = 300 K, 600 K, and 900 K. The graphs were created by merging snapshots at different time steps. Also shown are the interstitial atoms and their connections with the LSIAs at which they are located. The first point of observation is that the pair of interstitial atoms in LSIAs is in either BC or BS configuration for most of the time. (It is actually difficult to distinguish between BS and BC configurations for nonzero substrate temperatures.) However, the atom-LSIA-atom connection is not necessarily a straight line in the direction. The interstitial atoms exhibit superposition of thermal swinging and vibrational movements around their LSIA because of temperature. The atoms (not drawn for clarity) on the nearest lattice sites of the LSIA also thermally swing and vibrate around their lattice points. There are moments when the interstitial atoms and the atom on one of the LSIA candidates approximately line up. At these moments, an LSIA jump may occur, during which one of the interstitial atoms falls into the new LSIA and the other returns to the lattice point where the previous LSIA is located. After we quenched the simulation boxes to zero temperature, the SIA configuration that led to the migration of SIAs from one basal plane to another (denoted as PS’ in Ref. [17]) was indeed found. However, the appearance of this state is infrequent. Most of the time, the interstitial atoms tend to swing around their LSIAs at small angles. This leads to our second point of observation. In Fig. 2, it is seen that the LSIA exhibits one-dimensional migration at low temperature (T_s = 300 K). The LSIA ‘jumps’ forward and backward in a direction. With increasing temperature (Ts = 600 K), changes in the direction of LSIA movement are observed. LSIA may change its movement from one direction to another equivalent direction in the same basal plane, an off-line jump as defined in the previous subsection. It may also first jump to one of the two nearest 0001} planes and then quickly shiftto the direction (the off-plane jump). However, the probabilities for off-line jumps and off-plane jumps are significantly less than those of in-line jumps. The trajectories of LSIAs are thus observed to be constituted by line segments of various lengths. With the further increase of temperature (T_s = 900 K), the occurrence probabilities of off-line jumps and off-plane jumps increase, and thus the average length of the line segments decreases. Even so, the events in which an LSIA jumps forward or backward in a direction still occur most frequently.

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. (color online) Typical trajectories of LSIAs for temperature Ts = 300 K, 600 K, and 900 K (from top row to bottom row). Left column are views in [0001] direction, right column are views in direction. Dark red dot: hcp lattice point; circle: LSIA; gold dot and line: an interstitial atom and its connection with its LSIA; green dot and line: the same as gold dot and line but for another interstitial atom at the LSIA.

To conduct a more quantitative analysis of the migration features of Zr SIA, we calculated occurrence frequencies of LSIA jump events, , where the superscript α denotes the type of jump (i.e., in-line, off-line and off-plane), instead of calculating the diffusion coefficient of SIA. The jumping frequencies are more essential than the diffusion coefficient. The diffusion of continuous diffusion theory can be calculated from , where l_(α) is the jump length of type α. In addition, the occurrence frequencies can be used for time-step sampling in kinetic Monte Carlo simulations of defect evolutions.^[30]

We determined the occurrence of an LSIA jump event by the approach described in Subsection 2.2. If is the number of times event αoccurs in all 1000 simulation boxes in the 30-ps simulation time, then in unit s⁻¹ is calculated from . It should be mentioned that to correctly count the number of jump events, the time step should be small enough to avoid long distance jumps. For the recording time step (0.1 ps) used, only the jumps between two adjacent hcp lattice sites were observed.

Figure 3 shows as a function of the temperature. Also shown is the summation of the occurrence frequencies of all types of events . Because the total number of jump events is 19340 for T_s = 300 K and more for higher T_s. a statistical error of is estimated at less than 1%. From Fig. 3, it is seen that the in-line jump frequency is significantly higher than the off-line jump frequency and the off-plane jump frequency at most considered temperatures. At room temperature (T_s = 300 K), and are very small compared with . With increasing temperature in a range of 300 K ≤ T_s ≤ 600 K, increases but at a declining rate. For T_s ≥ 600 K, slightly decreases with increasing temperature. Conversely, and increase with rising temperature at an increasing rate in the range of 300 K ≤ T_s ≤ 600 K. For T_s ≥ 600 K, and tend to linearly increase with rising temperature and increases slightly faster than . The rates of the increase in and compensate for the rate of decreasing in . This leads to good linear dependence on temperature for the occurrence frequencies of all events. The dependence of on the temperature was very well fitted by the linear equation , with and . The linear temperature dependence suggests that the occurrence frequencies of all LSIA jump events follow the Einstein-Smoluchowski relationship and thus that LSIA migration tends to Brownian-like motion at our considered temperatures. However, the occurrence frequency for a specific type of event, for example, , may not always follow the Einstein–Smoluchowski relationship in the considered temperature range. Because the occurrence frequencies are different for the various types of LISA jumps, the Brownian-like motion of LSIA is anisotropic.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) Occurrence frequencies of event types, , as a function of temperature, where ■ denotes total frequency of events ; ● denotes ; ▲ denotes ; ◀ denotes ; and ▼ denotes . The red solid line is the linear fitting of : , with Ts in unit K and fE(Ts) in unit ps−1.

For comparing our results with other authors’, the occurrence frequency of in-plane jumps, , is also shown in Fig. 3. In contrast to , the dependence of the in-plane jump frequency on T_s should be considered in two temperature regimes divided by T_s = 600 K. This is in agreement with what was observed in Osetsky et al.’s work,^[10] but the in-plane jump frequency displayed in Ref. [10] is systematically larger than ours. Osetsky et al. have fitted the in-plane jump frequency to the Arrhenius relationship in two temperature regimes (T_s ≤ 600 K and T_s ≥ 800 K), with activation energies equal to 0.007 eV and 0.028 eV, respectively. The small activation energy in the low-temperature regime (0.007 eV), which is comparable to or even smaller than the considered temperature, is actually an indication that the SIA migration tends to be non-activated migration. We have implemented the Nudged–Elastic–Band (NEB) method^[31] to perform static calculation of diffusion barrier and the results are shown in Fig. 4. It is observed that the ILJ has the lowest diffusion barrier which is less than 0.01 eV (equivalent to about 100 K). This proves again that glide in direction should be the dominating mechanism of SIA migration at low temperature. The diffusion barriers of OLJ and OPJ are shown in Figs. (4b) and (4c). The OLJ path actually consists of an in-line glide and a rotation path with a diffusion barrier of 0.18 eV. The diffusion barrier of OPJ is 0.23 eV. The NEB results explain the dynamical result shown above.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) NEB analysis results of ILJ, OLJ, and OPJ. In ILJ, the initial state is set at L1 and the final state is set at L3. In OLJ, the initial state is set at L1 and the final state is set at L4. In OPJ, the initial state is set at L2 and the final state is set at L5, where the Ln is denoted in Fig. 1.

Figure 5 displays the percentages of occurrence frequencies of event types, defined by . Again, can be fitted by the linear equation . The fitting parameters A^(α) and B^(α) are given in Table 1. And in larger scale simulation, such as KMC, the occurrence frequency of type α can be calculated by .

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) Plots of percentage of occurrence frequencies of event types, , versus temperature. ● denotes ; ▲ denotes ; ◀ denotes . The solid line are the linear fittings of : , with the fitting parameters given in the following Table 1.

Table 1.

Table 1.

Table 1. Parameters for fitting to the percentage of occurrence frequencies of event types in Fig. 5. . ILJ OLJ IPJ OPJ A(α) –0.06477 0.04029 –0.02448 0.02448 B(α) 112.60166 –3.87286 108.7288 –8.7288

Table 1. Parameters for fitting to the percentage of occurrence frequencies of event types in Fig. 5. .

An LSIA has two nearest sites for ILJs and four nearest sites for OLJs (see Fig. 1). For the convenience of description, we call the migration fraction-dimensional if , , and . In contrast, we denote ‘full’ three-dimensional migration as the case in which the probabilities for an LSIA jumping to all of its 12 nearest neighbouring sites are the same, namely, satisfying ; ‘full’ two-dimensional migration as the case in which the probabilities of jumping to its 6 nearest neighbouring sites in the same basal plane are the same, namely, ; and ‘full’ one-dimensional migration as the case in which only in-line jumping exists, namely, and Clearly, the LISA migration in Zr displayed above is fraction-dimensional. In the next section, we will show that the kinetic effects of fraction-dimensional migration could be different from those that would occur if the migration of SIAs are considered to be full one-dimensional or full two-dimensional.

As demonstrated in the previous section, the migration of LSIA in the basal plane of Zr is fraction-dimensional. According to the random walk theory,^[32] a particle randomly walking on one-dimensional lattice sites can repeatedly visit the lattice sites many times. If off-line jumps occur, the visits to these sites are branched, and visits to sites aligned in other directions, which could not be visited when the particle walks in one direction, are probable. During their migration, SIAs can be annihilated when they encounter, for example, vacancies that are randomly distributed in the material. The SIAs clusters can form when the SIAs encounter each other. Obviously, the probability for an SIA encountering other defects is proportional to the number of different lattice sites that the SIA, or the equivalent LSIA, can visits. If N_V is the number of different sites that the LSIA visits after N_E jump events, we suggest using the ensemble-averaged value, , to measure the potential kinetic effects. It is difficult to obtain an analytical expression of n_SPE. We perform Monte Carlo (MC) simulations to calculate n_SPE. We start an MC simulation by randomly selecting an hcp lattice site for the initial LSIA, and mark this site as being visited. Then, the LSIA migrates through a series of jump events. For each LSIA jump, a random number rwith uniform distribution in the interval [0,1] is generated and the type of the jump event is determined according to the percentage of occurrence frequency . For example, if , the jump event is ILJ type; if the jump event is OLJ type; if , the jump event is OPJ type. With the determined event type α, the LSIA moves to the site randomly chosen among the candidate sites corresponding to the event type (see Fig. 1) The numbers of candidate sites are 2, 4, and 6 for ILJ, OLJ, and OPJ, respectively. A site can be visited for many times. However, a visited site can contribute only one count to the number N_V. The simulation terminates after issuing an assumed number N_E of events. N_E can be considered as the average number of jumps for a moving LSIA before its annihilation due to, for example, recombination between SIAs and vacancies. For a given N_E, many MC simulation runs are performed. n_SPE is then calculated by , with ⟨ ⟩ denoting the average over the runs.

Figure 6 displays two typical trajectories of LSIA migration obtained by the MC simulations with N_E = 10⁴, using , as given in Subsection 2.3, for T_s = 300 K and 600 K. The size of the simulation box in the MC simulations is 31 × 31 × 31 times the simulation box size used in the MD simulations above. Thus, the trajectories are much longer than what observed in MD simulations (see Fig. 2). The colour indicates the number of times that a lattice site is visited by the LSIA. It is observed that because of the small and at low temperature (T_s = 300 K), the LSIA migrates in a direction for a long line segment before changing direction and that many lattice sites are visited several times (Fig. (6a)). With increasing temperature (T_s = 600 K), more lattice sites are visited, and on average the number of times that a lattice site can be visited decreases (Fig. (6b)).

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) Typical trajectories of LSIAs viewed in [0001] direction obtained by MC simulations. (a) Ts = 300 K; (b) Ts = 600 K. The colour represents the number of times that a site is visited by the LSIA.

Figure 7 shows the plots of n_SPE versus N_E for different substrate temperatures. Also shown are the plots of for full three-, two-, and one-dimensional migrations of LSIA on hcp-lattices, where the superscript d denotes the dimensionality. Keeping in mind that mathematically, n_SPE = 1 if N_E = 1 in any case, all n_SPE values are seen to monotonically decrease with increasing N_E due to the spreading out of visited lattice sites. Lower migration dimensionality corresponds to lower levels of spread. Thus, and provide the minimum and maximum boundaries for n_SPE of the migration phenomenon in the hcp-lattice, respectively. The difference between the minimum and maximum is 60% at N_E = 20. The difference widens with increasing N_E. The large difference implies that n_SPE is very sensitive to the migration dimensionality. Returning to LSIA migration in Zr, it is seen that n_SPE varies from 38% to 79% at N_E = 20 and from 10% to 70% at N_E = 10⁴, when the temperature T_s increases from 300 K to 1200 K. Another observation is that the increase in causes n_SPE to asymptotically approach to a constant for large N_E. For full two-dimensional and low-dimensional migration, the slope of n_SPE as a function of N_E is large even for N_E = 1000. The large slope is probably an indication that the kinetics of the LSIAs can be dependent on the concentration of the sites with which the SIAs can interact. In such a case, the use of concentration-independent kinetic quantities, such as a diffusion coefficient that originates from the theory of continuous diffusion, should be revalidated. We will address this issue more generally in the future. Here, we address another question of the usage of diffusion coefficients when the LSIA migration is fraction-dimensional.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. (color online) Plots of average number of visited sites per event, nSPE, versus the number of jump events, NE. Also shown are nSPE for full 3D, 2D, and 1D migrations. The definitions of full 3D, 2D, and 1D migrations are found in the texts.

In Refs. [10] and [11], anisotropy of Zr SIA migration is obtained by differentiating the diffusion coefficients in the basal-plane and along the c-direction of hcp crystal structure. In kinetic Monte Carlo simulations or mean field rate theory, the diffusion coefficients are used as inputs. The diffusion of Zr SIA is usually treated aseither full one-dimensional or full two-dimensional migration in the basal plane.^[4,5,33] However, our present results reveal that the in-basal-plane migration of LSIA in Zr is fraction-dimensional. To study the consequence of such treatments, we perform further MC simulations for three cases. The first case is the ‘real’ fraction-dimensional case, in which all are those obtained by MD simulations in Subsection 2.2. For the other two cases, the same is also adopted. However, we set and to treat the in-basal-plane migration of LSIA as being fully one-dimensional and set and to treat the in-basal-plane migration of LSIA as being fully two-dimensional. The diffusion coefficient D⁽²⁾ in the basal-plane is then calculated by , and the three-dimensional diffusion coefficient D⁽³⁾ is calculated by , with (x, y, z) denoting the displacement vectors of LSIAs after N_E (≫ 1) jumps. The time spent in N_E jumps, t, can be calculated by .

Figure 8 shows plots of n_SPE versus diffusion coefficients and , obtained with N_E = 10⁴ for all three cases. The temperature T_s ranges from 300 K to 1200 K. Because an LSIA jump in the basal plane always contributes the same square displacement , where a₀ is the basal lattice length, to the mean square displacement (MSD) no matter whether the jump is in-line or off-line, the three cases have the same diffusion coefficients D⁽²⁾ (T_s) and D⁽³⁾ (T_s) at the same temperature, as denoted by the dash line in Fig. 8. However, the n_SPE of the ‘real’ case deviates significantly from the values of n_SPE of the assumed one- and two-dimensional cases. This result suggests that the conventional diffusion coefficients cannot reflect the real underlying kinetics of LSIA migration in Zr, in which the in-basal-plane migration of LSIAs is neither fully one-dimensional nor fully two-dimensional but rather fraction-dimensional. Thus, the models in downstream applications for modelling the microstructure evolution in Zr, such as KMC or mean field rate theory, should account for the effects of the fraction-dimensional migration of SIAs to make accurate predictions.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. (color online) Plots of nSPE versus diffusion coefficient D(2) (denoted by ◼) and D(3) (denoted by ●) for full two-dimensional case (red) and fraction-dimensional case (blue), and full one-dimensional case (black). See text for definitions of the cases.