Self-Diffusion in the Hexagonal Structure of Zirconium and Hafnium. Computer Simulation Studies

Instituto de Tecnología Prof. J. Sabato, Universidad Nacional de San Martín, Pcia. de Buenos Aires, Argentina Departamento Materiales, Centro Atómico Constituyentes, Comisión Nacional de Energía Atómica, Avda. Gral. Paz 1499, C.P. 1650, San Martín, Pcia. de Buenos Aires, Argentina, and Consejo Nacional de Investigaciones Científicas y Técnicas, Argentina Departamento Materiales, Centro Atómico Constituyentes, Comisión Nacional de Energía Atómica, Avda. Gral. Paz 1499, C.P. 1650, San Martín, Pcia. de Buenos Aires, Argentina


Introduction
A detailed revision on bulk and interface boundary diffusion in group IV hexagonal close packed (hcp) metals (Ti, Zr and Hf) and alloys has been recently carried out by Herzig et al. 1 .Also, an overview on self-and impurity diffusion in hcp Ti and Zr, due to Pérez et al., is found in quotation 2 .It is well known that metals of the IV group exhibit phase transition from α (hcp) to β body centered cubic (bcc), being the temperature of the reaction 1156 K for Ti, 1138 K for Zr and 2016 K for Hf.The relatively low transition temperature for Ti and Zr makes difficult diffusion experiments in the α-phase, while diffusion in α-Hf can be studied in a larger temperature interval.Considerable amount of experimental results collected in Herzig et al.
Herzig et al. 1 and Pérez et al.
Pérez et al. 2 support the fact that small metallic atoms, such as Fe, Cr or Ni, act as fast-diffusing impurities and produce a pronounced enhancement effect on both self-and substitutional bulk solute diffusion in these metals.Regarding self-diffusion behavior, the studies in α-Hf demonstrate generally a good linearity of the Arrhenius plot, and only at the lowest temperature a slight upward deviation can be supposed 1 .Also α-Ti, contrary to α-Zr, shows an Arrhenius behavior independently of the impurity content [1][2] .According to Herzig et al.
Herzig et al. 1 , self-diffusion in high purity materials is intrinsically 'normal' and is represented by linear Arrhenius plots which correlate well with selfdiffusivities in other hcp metals with no phase transition in the solid state.Precisely, α-Zr and α-Hf impurities free materials are the object of the present work, where computer simulation techniques are applied in order to analyze normal vacancy-mediated bulk diffusion.
The scheme of the present paper is as follows.In Section II, the calculation procedure is briefly described and the potential used to represent the interatomic interactions in α-Hf is presented.In Section III, the effect of the hcp lattice anisotropy on self-diffusion is considered.In Section IV, the obtained results for the vacancy migration are reported and the vacancy/tracer diffusivity in the hcp structure is discussed.Finally, in Section V, the calculated results are compared with the measured ones.

Calculation Procedure
Molecular static (MS) and molecular dynamic (MD) are used to calculate defect energies in crystallites where the atomic interactions are represented by many-body potentials of the embedded atom type (EAM) [3][4] .
MS simulations, based on the conjugate gradient method 5 , are carried out in crystallites of approximately N = 2000 atoms under constant volume condition.The quasi-static approximation 6 is applied to calculate the barrier for the vacancy migration.Briefly, a nearestneighbor divacancy conveniently oriented is created and an extra lattice atom is located in different positions between the vacant sites.In order to prevent recombination, no displacement towards the nearest vacancy is allowed.The (3N-1) coordinates that remain free are allowed to relax in order to minimize the energy under the interatomic potential, and in this way, to find the minimum energy path.
MD simulations are performed in order to study the temperature dependence of the tracer diffusion coefficient D t .Temperatures range from 800 to 1500 K for α-Hf and from 700 to 1000 K for α-Zr.In the absence of free surfaces or an appropriate vacancy density, the crystals remain stable within the times of the simulations even at the high temperature limits.A crystallite composed of N = 383 particles and a vacant site is simulated under periodic boundary conditions.For each temperature the system is balanced, as well as a reference perfect lattice, to null pressure in a canonical ensemble.Then, the evolution of the system at constant energy and volume is allowed for the period of time that assures a reasonable jump statistics.
The interatomic potential used to represent α-Zr has been developed by Pasianot and Monti 7 and previously employed in studies on the energetic of bulk vacancy and self-interstitials and the structure and energetic of grain boundaries 8 and surfaces 9 .The potential for α-Hf has been developed according to the same procedure followed for α-Zr.
The EAM scheme [3][4] assigns to each atomic site ' i ' an energy E i given by where V is a pair potential, R ij is the distance between atoms 'i' and 'j', F is the so-called embedding function and ρ i is an 'electronic-like' density given by a superposition of atomic pair-like functions In order to determine V, F (ρ) and ϕ, references 7 and 10 are followed.F (ρ) is imposed to have null first derivative at the perfect lattice density ρ 0 , this makes V (R) an effective pair potential.The function ϕ is chosen as a Thomas-Fermi like screening function, continuous up to the second derivative and smoothly matched to zero at the cut-off distance (see Equation 5and Table 1 in Pasianot and Monti 7 ).V (x) is a seven-piece cubic polynomial continuous up to the second derivative according to where H (x) is the Heaviside step function, x k stand for adequately chosen knot points and A k are coefficients determined in the following way.The seven coefficients plus the embedding function and its second derivative at the perfect lattice constitute nine unknowns exactly solved for through a system of equations involving the two lattice parameters a and c, the cohesive energy E c , the five elastic constants and an approximate value for the vacancy formation energy E f ν .The obtained coefficients for the α-Hf potential are reported in Table 1.
Notice that the potential for α-Hf reported in Pasianot and Savino 10 has a longer range that the present one, for which the x k 's are the same that those for the α-Zr potential.Finally, F (ρ) in Equation 1is obtained as indicated in Pasianot and Monti  Pasianot and Monti 7 .

Diffusivity in hcp Lattices
Diffusion by the vacancy mechanism in the hcp lattice takes place by two different first-neighbor jumps 6,11 .The tracer can either jumps to a vacant site located in the same basal plane, or to one located outside that plane.They are hereafter named b-jump and nb-jump, respectively.In terms of these jumps, the tracer diffusion components parallel and perpendicular to the c-axis are written as Above, C ν is the equilibrium vacancy concentration at temperature T given by where G f ν is the vacancy formation free energy and K B is the Boltzmann constant.Γ is the vacancy jump frequency for basal or non-basal jumps given by where G m ν is the vacancy migration free energy, and ν is the atomic attempt frequency.
Finally, f's stand for the correlation factors.For self-diffusion in the hcp structure, these factors are shown to depend on the ratio according to the following expressions 12 : where S m is the vacancy migration entropy, can also be evaluated for each type of jump.

Results
According to the quasi-static method applied in the MS simulations, the results for the vacancy migration energy quoted in Table 2 are obtained.
It is noted that, for both materials, basal jumps are slightly less energetic than non-basal jumps.
Regarding MD simulations, the atomic migration contribution to diffusion parallel to the c-axis is calculated as where t is the run length, < Z 2 (t) > is the mean-square displacement projected along the c-axis for all the atoms and B is the asymptotic value of the same quantity for a perfect bulk, which is independently calculated and subtracted in order to retain only the diffusive part 13 .Similar expressions are valid for directions X and Y, that are perpendicular to the c-axis.The associated migration energy is obtained from the slope of the D * (T) plot vs. 1/T.Particularly, the energy for the non-basal jump is obtained from the plot of D * II , while from the plots for directions X or Y an effective energy that corresponds both to basal and non-basal jumps is determined.In addition, the self-diffusion coefficients of Equation 4 are related to the corresponding D * (T) by where 1 / N is the vacancy concentration in our model system.
Figure 1 is the Arrhenius plot for D * II (T) calculated in both materials, the corresponding migration energies being E m nb = 0.57 ± 0.03 eV for α-Zr and E m nb = 0.79 ± 0.03 eV for α-Hf.Another way to determine the defect migration energy is from the jump frequency Γ.By MD, the jumps can be identified by checking, at each time step, the occupation number of elementary cells in a reference lattice.Every change in the occupation number implies that a jump has occurred.By recording the involved lattice sites and the occurrence time step, the jumps can be easily classified and the corresponding frequencies evaluated 13,14 .In doing that, the oscillations (they occur when an atom jumps and immediately returns to the original site) are filtered for the time interval τ ≤ 0.3 ps for α-Zr and τ ≤ 0.4 ps for α-Hf.They are approximately twice the corresponding Debye periods 13 .The oscillations would produce an increase of the jump frequency and no effect on the atomic diffusion.On the other hand, the return jumps that take place for a greater time interval than two Debye periods are included in Γ but discounted through the correlation factors.Figures 2a and 2b show the Arrhenius plot for the frequency of basal (solid line and black circles) and non-basal (dotted line and white circles) jumps for α-Zr and α-Hf, respectively.The slope of the straight lines gives the vacancy migration energy, the obtained values are quoted in Table 3.
Note that the obtained values for the non-basal jumps are in good agreement with those obtained from the mean square displacement projected along the c-axis and, as previously predicted by MS, basal jumps are also found to be slightly less energetic than non-basal ones.
Finally, the vacancy migration frequencies ν * defined in ( 8) can be obtained from the intercepts in Figures 1 and 2. Table 4 collects the results.
Values predicted by the mean square displacement are slightly higher than those obtained from the jump frequency, and non basal jumps have higher vacancy migration frequency than the basal ones.Taking ν in (8) as the Debye value, it is seen from Table 4 that the vacancy migration entropy results in the range 0.9 to 1.6 K B , in reasonable agreement with the values deduced from experiments at least for fcc metals 15 .
The results obtained in this work allow the evaluation of the diffusion activation energy Q = E f ν + E m ν for vacancy-mediated bulk mechanism.The vacancy formation energy E f ν , which unrelaxed value is a parameter fitted by the potential, is predicted to be 1.74 eV for α-Zr and 2.11 eV for α-Hf.
As previously discussed in Pasianot and Monti 7 , the above E f ν value for α-Zr could be lower than the experimental one.According   to Hood et al. 16 , the positron annihilation spectroscopy technique (PAS) gives only a lower bound of about 1.5 eV, so a higher value than the one fitted by the potential could be expected.In this line, by applying the full-potential linear muffin-tin orbital method, an unrelaxed E f ν value of 2.07 eV 18 is predicted.On the other hand, the present results for E m ν , although slightly dependent on the simulation method, are consistent with experimental findings based on PAS that estimate a value in the range 0.6-0.7 eV 17 .
Finally, the diffusion activation energy here calculated is about 2.30 eV 7 .This result is lower than the ones obtained from experiments.Effectively, 3.06 eV is reported in Herzig et al. 1 and attributed to an intrinsically normal behavior, while in Pérez et al. 2 , and for temperatures below 900K, 3.50 eV is reported from a non-linear Arrhenius plot.Assuming a normal diffusion behavior, we consider that the discrepancy between experiments and calculations can be a consequence of the E f ν value fitted by the potential, as previously discussed in Pasianot and Monti 7 .
In α-Hf, and following a common practice 10,19 , the interatomic potential has been fitted to an unrelaxed vacancy formation energy taken as a third of the cohesive energy.That value, 2.15 eV, is in reasonable agreement with the one obtained from ab initio calculations, 2.37 eV 18 , and with the experimental result, 2.45 ± 0.2 eV 20 .
No measured values for the vacancy migration energy are known to the authors, the ones calculated in this work being around 0.8 eV.Finally, the calculated diffusion activation energy is about 2.93 eV.This prediction seem to be in reasonable agreement with the experimental result of 3.23 eV 1 .

Conclusions
It has been shown a good agreement between results obtained by Molecular Static and Molecular Dynamic in relation to the vacancy migration parameters in α-Zr and α-Hf modeled by EAM type interatomic potentials.
For both materials and computer simulation techniques, the vacancy migration energies and entropies for basal jumps are slightly lower than those for the non-basal ones, in agreement with the early Zener's predictions 21 that larger migration energy implies also larger migration entropy.
For α-Hf, the interatomic potential presented in this work predicts a diffusion activation energy in reasonable agreement with the value obtained from experiments in Herzig et al. 1 .

Figure 1 .
Figure 1.Atomic migration contribution to the diffusion coefficient D * II obtained from Equation9for lattices with one vacancy: α-Zr (black circles) and α-Hf (white circles).

Figure 2 .
Figure 2. Temperature dependence of the jump frequencies.
The purpose of this work is to determine the vacancy migration energies associated to both b and nb-jumps, indicated as E m b and E m nb , respectively, and compare the values predicted by MS and MD simulations.The quantity known as vacancy migration frequency

Table 1 .
V (x) (in eV) potential parameters, x k in units of a.

Table 2 .
Vacancy migration energy calculated by molecular static.

Table 3 .
Vacancy migration energy calculated from the jump frequency evaluated by molecular dynamic.

Table 4 .
Vacancy migration frequency obtained from the Arrhenius plots of both the mean square displacement (in parenthesis) and the jump frequency.