University of Birmingham Diffusion anisotropy of poor metal solute atoms in hcp-Ti

Atom migration mechanisms influence a wide range of phenomena: solidification kinetics, phase equilibria, oxidation kinetics, precipitation of phases, and high-temperature deformation. In partic-ular, solute di ff usion mechanisms in α -Ti alloys can help explain their excellent high-temperature behaviour. The purpose of this work is to study self- and solute di ff usion in hexagonal close-packed (hcp)-Ti, and its anisotropy, from first-principles using the 8-frequency model. The calculated di ff usion coe ffi cients show that di ff usion energy barriers depend more on bonding characteristics of the solute rather than the size misfit with the host, while the extreme di ff usion anisotropy of some solute elements in hcp-Ti is a result of the bond angle distortion. C 2015 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 Unported License. [http: // dx.doi.org / 10.1063 / 1.4921780]


I. INTRODUCTION
The study of diffusion mechanisms is fundamental to understanding how materials behave in a variety of conditions. The migration of solute atoms, such as Si, is particularly important to explain the high-temperature behaviour of α-Ti and near-α-Ti alloys since this behaviour is largely controlled by dislocation climb. Similarly, the migration of solute atoms also controls ageing and precipitation kinetics, both during processing and while in service. The study of Si diffusion mechanism is important for high-temperature Ti alloy design. This study compares the diffusion behaviour of Si, usually added to near-α-Ti alloys to improve their creep properties, to elements located close to it in the periodic table (Al, Ga, Ge, In, and Sn). Such an approach can highlight trends and help explain relevant mechanisms.
At low temperatures, below 882 • C, Ti maintains a hexagonal close-packed (hcp) structure (α-Ti). As a result of the hcp crystal structure, diffusion can be highly anisotropic, with diffusivity along the basal plane (denoted in this article as D ∥ ) being faster than diffusivity along the normal direction to the basal plane (denoted in this article as D ⊥ ). Therefore, it is not trivial to study diffusion mechanism in α-Ti experimentally. These experiments would require the use of single-crystal α-Ti samples, which is hard to achieve. Indeed, most experimental measurements of diffusivity in α-Ti focus on determining a bulk diffusion coefficient in poly-crystalline samples. Experimental data on the diffusion of poor metals (Si, Al, Ga, In, Sn) in α-Ti are scarce and they are summarised in Table I. Experimental studies of Ge diffusion in hcp-Ti could not been found in the literature.
Köppers et al. have studied Al, 1,2 Ga, In, 2 and selfdiffusion 1 in ultra-pure single crystal α-Ti. They have measured the diffusivity within the basal plane and fitted the values for pre-exponential factor, D 0∥ , and diffusion activation energy, Q ∥ . The data of Al-and self-diffusion coefficients a) Electronic mail: lxs234@bham.ac.uk b) Electronic mail: a.mottura@bham.ac.uk perpendicular to the basal plane have allowed them to calculate an anisotropy factor equal to 0.63 and 0.5, respectively. Other studies on the diffusion of Sn 3 and Si 4 in α-Ti have been conducted without taking anisotropy into consideration. Therefore, a detailed understanding of the diffusion mechanisms of Si, Ga, Ge, In, and Sn in α-Ti is still not present.
In the last few years, first-principles calculations have been successfully used to evaluate self-diffusion 5 and solutediffusion 6 in hcp-Mg and hcp-Zn, and O diffusion in hcp-Ti. 7 These show that density functional theory (DFT) is a valid instrument to determine diffusion parameters, as well as helping our understanding of complicated diffusion mechanisms. The purpose of this work is to extend these firstprinciples studies to the poor metals diffusion in α-Ti. This will be conducted using DFT as implemented in the Vienna ab initio simulation package (VASP) 8 and the climb image nudged elastic band (CI-NEB) method. 9

II. METHODS
All calculations were obtained using DFT as implemented in VASP 5.3.2. 8 The projector augmented wave (PAW) 10 method was used to describe the electronic structure and the exchange and correlation functional was computed using the generalised gradient approximation (GGA) parametrised by Perdew-Burke-Ernzerhof (PBE). 11 The accuracy was guaranteed by a Monkhorst-Pack mesh of 32 × 32 × 17 k-points for a 2-atom unit cell (scaled to maintain k-point density for larger supercells) and an energy cutoff of 350 eV. Defect energies such as vacancy formation energies, E v , and binding energies of solute-vacancy pairs, ∆E b , were calculated using 3 × 3 × 3 supercells (54 atoms).
The choice of this size of supercell was based on the best balance between accuracy and computational time. Vacancy formation energies calculated with 3 × 3 × 3 and 4 × 4 × 4 supercells (128 atoms) are similar, respectively, 1.976 eV and 1.971 eV, while it is equal to 1.964 eV in case of a smaller supercell of 2 × 2 × 2 (16 atoms). I. Diffusion parameters of Al, Ga, In, Si, Sn diffusion and self-diffusion in α-Ti. Q ∥ and D 0∥ are the activation energy and pre-exponential factor of diffusion within the basal plane, respectively; D ⊥ /D ∥ is the diffusion anisotropy ratio; Q and D 0 are the activation energy and pre-exponential factor of bulk diffusion. The GGA for the exchange and correlation functional was chosen because it is more accurate for evaluating bulk properties such as lattice constants, elastic constants, and cohesive energies 12,13 when compared to the local density approximation (LDA). 14 Both exchange and correlation functionals underestimate the vacancy formation energy relative to experiments. 15 The error is usually higher in case of the GGA; however, both approximations to the exchange and correlation functional give similar results in the case of hcp-Ti. We find that E v is equal to 1.985 eV using the LDA and 3 × 3 × 3 supercell, similar to the results obtained with the GGA. This is also confirmed by E v values for hcp-Ti evaluated by Medasani et al. 16 using both the LDA and GGA.
The NEB method 9 allows to determine the minimum energy pathways minimising and interpolating the energy of a set of images between the initial and final configurations. In this work, the NEB was implemented with climb image method (CI-NEB) 9 and improved tangent estimation. 17 The CI-NEB calculations were performed using same setup (supercell size, k-points mesh, and energy cutoff) used for defect energy calculations.
Diffusion coefficients of hcp crystals can be evaluated using the 8-frequency model developed by Ghate in 1964. 18 In his theory, Ghate defines the diffusivity, D, as a function of frequency jumps. As discussed previously, the hcp anisotropy leads to the diffusion within the basal plane being different from the one perpendicular to the basal plane. For simplicity, we use the notation ∥ for diffusion within the basal plane and ⊥ for diffusion perpendicular to the basal plane. These two jumps are schematically represented in Fig. 1. D ∥ and D ⊥ are defined as 18 where a and c are lattice parameters, C is the probability of having a vacancy sitting beside a solute atom, f ∥x , f ⊥x and f ⊥z are partial correlation factors, and w ∥ and w ⊥ are impurity jump frequencies schematically represented in Fig. 1.
The probability of finding a vacancy sitting next to a solute atom is calculated using the formation energy of a vacancy next to the solute, ∆H v , E v depends on the energy of perfect cell and a cell containing a vacancy while ∆E b is defined as the energy difference between a system containing the vacancy-solute pair and a system where the solute and vacancy are at infinity In Eqs. (4) and (5), M and X are notations for the solvent and solute atom, respectively, and N is the total number of atoms. Therefore, pair, respectively. Note that the solute-vacancy binding energy is positive when repulsive and negative when attractive.
The correlation factors, f ∥x , f ⊥x , and f ⊥z , depend on two impurity-jump frequencies (w ∥ , w ⊥ ) and six solvent-jump , as outlined in Fig. 1, where λ ∥ is the jump distance within the basal plane, λ ⊥B is the projection of the w ⊥ jump distance on the basal plane. S ∥x and S ⊥x are calculated solving the following equations: where F is chosen to be equal to 0.736 18 and quantifies the effect of vacancies returning back to the nearest-neighbour positions. Each w is a function of its migration energy, ∆H m , and effective frequency ν * , 19 according to assuming Vineyard's definition of ν * , 20 where ν i and ν ′ j are normal frequencies of ith degree of freedom of initial configuration and jth degree of freedom of saddlepoint configuration. Only frequencies of the jumping atom were considered. The migration energy, ∆H m , is defined as the difference between the energy at the transition state, E TS , and the energy at the initial state, E IS , along the minimum energy path.
In the case of self-diffusion, w a , w ′ a , w c , and w ′ c are equal to w ⊥ , while w b and w ′ b are equal to w ∥ . Therefore, only these two frequencies are required to calculate the self-diffusivity.

III. RESULTS
The vacancy formation energy in pure α-Ti is equal to 1.976 eV and is in good agreement with previous firstprinciples studies. 19,[21][22][23][24] The binding energies of vacancysolute pairs for Al, Si, Ga, Ge, In, and Sn in hcp-Ti are shown in Fig. 2. The binding energy when vacancy and solute are in adjacent basal planes, ∆E b ⊥ , is displayed using black bars. The grey bars indicate the binding energy when vacancy and solute are in the same basal plane, ∆E b ∥ . A positive value indicates a repulsive force between vacancy and solute. On the contrary, ∆E b is negative when there is an attraction between vacancy and solute. There is an energy benefit in having a vacancy located adjacent all elements except for Al and Ga.

A. Self-diffusion
The migration energies of self-diffusion along ⊥ and ∥ directions were evaluated using 3 × 3 × 3 supercells and 5 images with CI-NEB. In case of the ∥ jump, ∆H m is equal to 0.413 eV, while ∆H m⊥ is equal to 0.427 eV. The perpendicular jump is less energetically favourable. Activation energies differ by only 6 meV: Q ∥ is 2.393 eV and Q ⊥ is 2.399 eV. Selfdiffusivity (D ∥ and D ⊥ ) in the temperature range of 873-1133 K is shown in Fig. 3. Ti diffuses more slowly normal to the basal plane than within the basal plane with anisotropy ratio D ⊥ /D ∥ of 0.89. D 0∥ is 0.978 × 10 −6 m 2 /s, while D 0⊥ is equal to 0.934 × 10 −6 m 2 /s. The present DFT calculation of self-diffusivity of hcp-Ti is compared to the first-principles calculations of Shang et al. 22 and experimental survey of Köppers et al. 1 Fig. 3. Shang and co-workers have used the GGA-PBEsol approximation and PAW method. Their calculations were done using a Monkhorst-Pack mesh of 11 × 11 × 11 k-points and 500 eV energy cutoff with a 3 × 3 × 2 (36 atoms) supercell. Köppers measured self-diffusion in ultra-pure Ti single crystals using radio-tracer 44 Ti.

in
Self-diffusion coefficients of the present work are consistent with the previous DFT values. However, Shang has calculated larger defect and activation energies (E v , Q ⊥ , and Q ∥ are equal to 2.14 eV, 2.65 eV, and 2.57 eV, respectively), and preexponential factors are one order of magnitude higher (D 0⊥ and D 0∥ are 1.20 × 10 −5 and 1.07 × 10 −5 m 2 /s). The difference between vacancy formation energy values is probably due to the choice of approximation for the exchange and correlation functional. In fact, Medasani et al. 16 have shown that E v is 2.18 eV when calculated using GGA-PBEsol. All these approximations to the exchange and correlation functional appear to overestimate the experimental data for the vacancy formation energy. 25 The discrepancy could also be attributed to the different supercell size, k-point mesh, and energy cutoff. In this case, however, we found that ∆H m values increase only 0.05-0.06 eV when using 4 × 4 × 4 supercells. Similarly, increasing k-point density and energy cutoff does not affect the migration energies significantly: the values obtained using an energy cutoff equal to 500 eV and a higher density k-point mesh differ by 1 meV with respect to the values obtained when using the chosen k-point density and energy cutoff for this work.
Computational results overestimate the diffusivity when compared to the experimental data obtained by Köppers et al. 1 who have fitted D 0∥ and Q ∥ , finding them equal to 1.35 × 10 −3 m 2 /s and 3.14 eV, respectively. The limited amount of data for D ⊥ measured by Köppers et al. has not allowed the authors to calculate D 0⊥ and Q ⊥ , but allowed for an estimate for D ⊥ /D ∥ which is given as 0.5. Our first-principles results give an anisotropy ratio of 0.89. This agrees with the common trend, D ⊥ /D ∥ < 1, for the hcp transition metals.

B. Solute diffusion
The 8-frequency model 18 was used to evaluate the diffusion (D ∥ and D ⊥ ) of Al (see Fig. 4). Also in the case of Al, diffusion within the basal plane is faster with respect to the diffusion normal to the basal plane. The anisotropy ratio calculated is equal to 0.75, similar to the self-diffusion value. A comparison of the data obtained in this work with computational data obtained by Mantina 26  of Al diffusion is equivalent to 0.63. As shown in Fig. 4, our results are more consistent with the experimental study than with the values obtained by Mantina. The disagreements with other DFT results can be attributed to the fact that CI-NEB is more accurate than NEB and to the use of different approximations for the exchange and correlation functional, k-point mesh, energy cutoff, and supercell size. Finally, the 8-frequency model 18 was applied to the diffusion of Si, Ga, Ge, In, and Sn in α-Ti. The migration energies were evaluated using 3 × 3 × 3 supercells and 5 images with CI-NEB and compared with self-diffusion values (Fig. 5). The values for ∆H v , ∆H m and Q for each impurity are reported in Table II. Si has the highest activation energy, while In and Sn are energetically the most favourable diffusers. The jump within the basal plane is energetically favourable with respect to the jump normal to the basal plane for all atoms. D ⊥ and D ∥ for Al, Si, Ga, Ge, In, and Sn are compared with selfdiffusion in Ti in Fig. 6. Si is the slowest diffusing element, TABLE II. First-principles evaluation of the formation energy of vacancy adjacent to the impurity, ∆H v , migration energy, ∆H m , and activation energy Q of poor metal diffusion in α-Ti along ⊥ and ∥ direction.  while In diffuses faster than all other elements considered in this study. Table III summarises the exponential pre-factor values for each element. The anisotropy ratios for In, Sn, and Ga are 0.53, 0.31, and 0.22, respectively, smaller than self-and Al anisotropy ratio. In the case of Si and Ge, D ∥ is two-three order of magnitude higher than D ⊥ leading to very low D ⊥ /D ∥ values, equal to 0.053 and 0.008, respectively.

A. Atomic size and solute diffusion
Similar to the transition metals, 27 atomic size plays a minor role in the vacancy-mediated diffusion of poor metals in α-Ti. Comparing the diffusivities obtained, the largest atom (In, with atomic radius of 1.54 Å) is the fastest diffuser amongst the poor metals, while Si is the slowest diffusing poor metal solute in Ti, although it has a relatively small atomic radius (1.11 Å). In Fig. 7, the migration energy values are compared with the atomic radius; Si has the highest ∆H m and In has the lowest ∆H m . This behaviour can be explained analysing the electron charge density difference between crystals where Si and In are located next to the vacancy and crystals with a vacancy (Fig. 8). A charge accumulation can be observed around Si atom suggesting a strong solute-solvent bond. As a result, the energy barrier to move a Si atom to the vacancy site is higher. On the contrary, In is not surrounded by a charge density accumulation; hence, the In-Ti bond is weaker requiring less energy for the In atom to exchange with the vacancy only.

B. Diffusion anisotropy
Anisotropy ratios of poor metal diffusion in α-Ti are all less than one (Fig. 9). This means that diffusion within the basal  plane is faster than diffusion normal to the basal plane. Al has the highest anisotropy ratio, close to self-diffusion value; In, Sn, Ga, Si, and Ge have significantly smaller anisotropy ratio. Ge stands out, since its diffusivity within the basal plane is three orders of magnitude faster than diffusivity normal to the basal plane. This effect is primarily due to ∆H Ge m∥ being half of ∆H Ge m⊥ . In 1957, for the fcc crystal, Swalin 28 proposed that ∆H m consists of two terms: (a) the energy due to the host lattice dilation at the transition state, (b) the compressibility of the solute necessary to squeeze thorough the window formed by the nearest neighbours at the transition state. The windows formed by the host atoms at the saddle point in the case of the hcp crystal are shown in Fig. 10. The first term depends on the displacement of this window: a large displacement leads to a small energy and vice-versa. We can approximate this by considering the distance between the solute atom at the saddle point and the nearest-neighbouring relaxed host atoms (see red dashed circle in Fig. 10). The second term is related to compressibility of the solute itself and can be inferred by computing the transition state energy while keeping the host lattice unrelaxed. In the case of the hcp lattice, these are not the only two factors that influence ∆H m . Considering the case of Ge, dilation does not significantly differ between the two jumps: R ∥ is equal to 2.67 Å while R ⊥ is 2.69 Å. The difference in host lattice dilation between the within-plane and out-of-plane cases is small and comparable to the host lattice dilation observed for the other solute atoms considered in this study. Similarly, the compressibility of the solute atoms at the saddle point does not change significantly: ∆H m∥ is 2.444 eV while ∆H m⊥ is 2.269 eV.
Therefore, in the case of the hcp lattice, lattice dilation and solute compressibility are not the only two factors that influence ∆H m . A lower energy saddle point can be accommodated by a severe distortion of the bond angles at the saddle point (γ in Fig. 10). In the unrelaxed lattice, these angles are equal to 68.8 • and 71.1 • , for ∥ and ⊥ jump, respectively. Allowing relaxation, bond angles are 62.2 • for w ∥ and 66.6 • for w ⊥ . Keeping the ideal bond angles but letting the host window dilate leads to ∆H m∥ increasing to 1.121 eV, bringing it closer to ∆H m⊥ . In contrast, ∆H m⊥ increases from 1.018 eV to 1.068 eV. This confirms that bond angle distortion leads to a low-energy saddle point configuration in the case of the ∥ jump.

V. CONCLUSION
In conclusion, we have studied self-and solute diffusion in hcp-Ti from first-principles using DFT and the CI-NEB method. Self-diffusion and Al diffusion values are consistent with the previous DFT and experimental data.
Si is the slowest diffuser amongst the poor metals, while the largest atom, In, is the fastest. This is attributed to the higher activation energy for diffusion in the case of Si, which is due to charge accumulation around Si. The findings are in agreement with previous work on transition metal solute diffusion in fcc crystals. The low diffusion of Si may play an important role on dislocation climb and consequentially on the high-temperature properties of Ti. A detailed understanding of the dislocation climb mechanism would be required to put this on a firm quantitative basis.
The diffusion barrier anisotropy of hcp-Ti is a result of host lattice dilation and saddle point nearest neighbour configurations. These lead to lower energy configurations at the saddle point in the case of transitions within the basal plane. In the case of Ge, this leads to an extremely low ∆H m∥ /∆H m⊥ ratio.