Interstitial diffusion of O , N , and C in α-Ti from first-principles : Analytical model and kinetic Monte Carlo simulations

• Users may freely distribute the URL that is used to identify this publication. • Users may download and/or print one copy of the publication from the University of Birmingham research portal for the purpose of private study or non-commercial research. • User may use extracts from the document in line with the concept of ‘fair dealing’ under the Copyright, Designs and Patents Act 1988 (?) • Users may not further distribute the material nor use it for the purposes of commercial gain.


I. INTRODUCTION
The combination of low density and excellent creep behaviour of Ti alloys makes these materials desirable for aerospace applications, such as gas turbine engines.However, the low resistance against hot gas limits their upper usage temperature. 1This detrimental phenomenon is due to the high affinity, solubility, and diffusivity of light elements, O, N, and C, in α-Ti.These elements are "α stabilisers" (i.e., they promote formation α-Ti at high temperature) and they have very high solubility values in α-Ti: 30 at.%, 19 at.% and 2 at.3][4][5] In addition, these elements have high diffusion coefficients, 10 3 -10 8 times higher than selfdiffusivities, as shown in Fig. 1.At high temperature, they dissolve and diffuse very fast in the Ti matrix, causing the formation of very hard and brittle surface layer of α.The high diffusivity of O, N, and C has a significant effect on the growth of this phase, since the depth of alpha phase is proportional to the square root of the exposure time. 6Understanding the migration mechanisms of these elements in α-Ti is extremely important to prevent such dangerous phenomena and design high-temperature Ti alloys.
The anisotropic character of the hcp crystal requires to perform single-crystal diffusion studies in order to fully understand the mechanism.However, obtaining a singlecrystal sample is not easy.In addition, the experimental approach would not give information on fundamental aspects a) Electronic mail: LXS234@bham.ac.uk b) Electronic mail: a.mottura@bham.ac.uk of the atom migration, the transition pathways.Knowing which interstitial sites O, N, and C prefer to occupy and the paths through which they migrate is a fundamental step for understanding the α phase phenomena.
In the last few decades, there has been an increasing utilisation of first-principles calculations on the study of vacancy-mediated diffusion, [7][8][9][10] as well as interstitial diffusion 11,12 in hcp crystals.In Wu and Trinkle's and O'Hara's articles, 11,12 ab initio calculations were used to provide data for an analytical approach called the Multi-State Diffusion (MSD) method, developed by Landman and Shlesinger, 13,14 in order to evaluate the diffusivity of interstitial elements.The results published in these articles show the reliability of Density Functional Theory (DFT) and the multi-state diffusion method.
This paper investigates the migration mechanism of O, N, and C in α-Ti using DFT together with the MSD method and kinetic Monte Carlo (KMC).
where k b and T are the Bolzmann's constant and absolute temperature, respectively.∆H m is the migration energy, defined as energy difference between the transition and initial states.The minimum energy paths were determined by climb image-nudged elastic band method (CI-NEB) 23 and improved tangent estimation 24 using 5 images between initial and final states with relaxation of ionic positions as well as shape and volume of the cells.The effective frequency ν * was approximated using Vineyard's definition. 25

A. Multi-state diffusion method
The multi-state diffusion method is an analytical method describing the continuous-time random walk of an atom migrating through internal states 13,14 and it can be applied to study interstitial diffusion in crystals.It consists of defining the stable interstitial sites occupied by the solute in the solvent crystal, followed by mapping the n possible jumps for each stable interstitial site.Each jump is defined by the average displacement and "waiting-time distribution" function, which describes the temporal dependence of each jump.The following step is building the Laplace transform matrix of waiting-time function, − → ψ (u), and Fourier transform matrix of displacement, p * ( − → k ) .The diffusion variance, σ 2 r (t), along a general direction r can be obtained solving where ∆ is the determinant of the matrix, where I is the identity matrix, and ∆ 0 is defined as Finally, the diffusion coefficient along r is defined as, The full matrices are provided within the supplementary material at Ref. 26.

B. Kinetic Monte Carlo
Kinetic Monte Carlo 27 is a particular Monte Carlo 28 method used for processes with known rates such as atom migration.It consists of mapping N possible events that can occur from a given state.Each event is defined by jump frequency, w i , displacement, and cumulative function of the FIG. 2. Schematic representation of interstitial sites in the hcp crystal: octahedral (OC), tetrahedral (TE), crowdion (CR), split dumbbells along c axis (SP), basal octahedral (BO), hexahedral (HE), basal crowdion (BC), and split dumbbells in the basal plane (BS).
Reuse of AIP Publishing content is subject to the terms: https://publishing.aip.org/authors/rights-and-permissions.Downloaded to IP: 147.188.224.127On: Mon, 11 Jul 2016 09:51:34 jump frequency, R i .The main algorithm can be summarised in the following steps: 1. set initial time t = 0 and the initial state of the system; 2. list all N possible events and their rates w i ; 3. for each event, calculate the cumulative function R i , defined as 4. pick a random number u ∈ (0, 1]; 5. find and perform the event i that satisfies R j < uR N < R i ; 6. pick a second random number u ′ ∈ (0, 1]; 7. calculate 8. repeat from step 2 with the system in the new state i and The diffusivity is defined as, where d is the dimension of displacement vector, ⟨t⟩ and ⟨r 2 ⟩ are the average time and square displacement between the initial and final positions of the species considered of n simulations.In this work, the diffusion coefficients were obtained as an average value of 100 KMC simulations of 10 6 solute jumps for each temperature.

C. Interstitial sites in hcp
There are 8 interstitial sites in hcp-crystals: octahedral (OC), tetrahedral (TE), crowdion (CR), split dumbbells along c axis (SP), basal octahedral (BO), hexahedral (HE), basal crowdion (BC), and split dumbbells in the basal plane (BS), schematically displayed in Fig. 2. The interstitial energies were calculated using, where N is the total number of solvent atoms, E{M N + I} is the energy of system where solute occupies an interstitial site, E{M N −1 + S} is the energy of system with substitutional solute and E{M N } is the energy of the perfect site.This energy is the energy difference between the interstitial and substitutional point defect of a given atom.A negative value of E i indicates that the solute tends to occupy the interstitial site considered rather than sitting substitutionally within the hcp lattice.

III. RESULTS
The interstitial energies of O, N, and C in α-Ti were evaluated using DFT calculations.They are presented in Table I and displayed in Fig. 3.When the interstitial site is unstable and relaxes to a different site, the final configuration is reported in brackets in Table I.The pathway networks were mapped and calculated using CI-NEB as described in Section III A.

A. Interstitial sites in α-Ti
Oxygen stable interstitial sites are OC, HE, and CR as shown in Table I.TE, SP, BC, and BS, instead, are unstable and relax to the HE site.The basal octahedral is also unstable since it has one imaginary vibrational frequency.In the hcp unit cell, there are two equivalent OC sites, six equivalent CR sites, and two equivalent HE sites (Fig. 3 II summarises the migration energies and attempt frequencies of these jumps. Similar to oxygen, nitrogen is unstable at TE, SP, and BS sites and relaxes to the HE or BC sites.BO has one imaginary frequency, hence it is also unstable.The stable sites are OC, HE, CR, and BC, as displayed in Fig. 3(b).The migration pathway network of N contains all the jumps described for O plus the transition from and to the BC sites.There are six neighbour BC positions around OC sites, leading to as many OC → BC jumps.Around CR, there are four BC sites, two on the basal plane above and two on the basal plane below, leading to four CR → BC transitions.Finally, the BC sites sit next to four CR and two OC sites: four BC → CR, and two BC → OC jumps.The BC ↔ HE jumps occur through crowdion sites.
The migration pathways network of N is illustrated in Fig. 4 and Table III presents the migration energies and attempt frequencies.
Carbon is stable at OC, CR, and BC sites in α-Ti (Fig. 3(c)).HE and BO are unstable since they have an imaginary vibrational frequency.The other interstitial sites relax to other configurations (see Table I).The transition network is built with the same procedure used for O and N and are shown in Figs.4(a

B. Diffusivities of O, N, and C in α-Ti
The jump frequencies and transition networks were used to solve the MSD equation and find the diffusion coefficients.KMC simulations were also performed and used to calculate the bulk diffusion.Perpendicular and parallel diffusivities of O calculated with MSD and KMC are shown in Fig. 5.They are compared with previous ab initio study of Wu and Trinkle. 11N and C diffusivities are displayed in Figs. 6 and 7, respectively.Tables V and VI contain the estimated activation energies and pre-factor diffusion coefficients for both the approaches, MSD and KMC.The diffusion within the basal plane is faster with respect to perpendicular diffusion for all atoms considered.
The results of this work for oxygen are compared to the findings of previous studies.The diffusion coefficients of O have been evaluated computationally by Wu and Trinkle 11 (Fig. 5) and experimentally by Bregolin et al. 15 (Fig. 8).Wu and Trinkle have evaluated interstitial sites and jump rates using VASP with ultrasoft Vanderbilt type pseudopotentials approximation 29,30 and generalised gradient approximation of Perdew and Wang 31 with ten valence electrons, energy cutoff of 400 eV, and k-point mesh of 2 × 2 × 2 for 96-atom supercell and using the CI-NEB method to estimate the transition state. 11hey have evaluated Q ⊥ and Q ∥ to be equal to 2.084 eV and D 0⊥ and D 0 ∥ to be equal to 2.18 × 10 −8 m 2 /s.Bregolin The main differences between our results and the results of Wu and Trinkle are that Wu and Trinkle performed CI-NEB calculations with a constant-shape supercells, describing the electronic structure with PAW potentials that treat 3p orbitals as valence rather than core electrons, and using a higher energy cutoff value.In order to investigate the influence of these parameters, the migration energies were calculated using the setup described in Section II with constant-shape supercells, ∆H cs m , as well as using the same energy cutoff (400 eV) and supercell size (96-atom) utilised by Wu and Trinkle allowing full relaxation, using the four-valenceelectron and ten-valence electron potentials, ∆H 400-96 m and ∆H 400-96-pv m , respectively.These results are compared with the calculations of Wu and Trinkle in Table VII.The choice of higher energy cutoff and larger supercell appears to have a minor impact, since ∆H 400-96 m values are similar to the migration energies obtained with smaller cutoff energy and supercell (see Table II).The migration energies calculated including the additional valence electrons (∆H 400-96-pv m ) are closer to the results of Wu and Trinkle when compared to the results obtained without treating the 3p electrons as valence electrons (∆H 400-96 m ).Nevertheless, the migration energies obtained keeping the supercell shape fixed (∆H cs m ) replicate the work of Wu and Trinkle well, especially in case of CR → OC and CR → HE jumps.Additionally, the choice of PAW potential has negligible effects on the defect energies of the HE and CR sites.The energy differences between the OC site and the HE and CR sites calculated with the 4s PAW are, respectively, 1.24 eV and 1.40 eV, similar to the results obtained with the 3p PAW calculations, 1.18 eV and 1.57 eV, respectively.This indicates that the discrepancies between the values of this work and the results of Wu and Trinkle are mainly due to the fact that our calculations have been performed with a full relaxation of the supercells.Due to this discrepancy, Wu and Trinkle have opted to assume that CR sites thermalise at high temperature, resulting in OC → OC jumps within the basal plane and HE → HE jumps, and a simplified equation for D ⊥ and D ∥ .Based on the transition pathways, this approximation mostly affects the evaluation of parallel diffusion.We have tested this by replacing OC → CR → OC paths with OC → OC paths and HE → CR → HE paths with HE → HE paths in the KMC simulations, obtaining diffusivities closer to those using Wu and Trinkle's formula if inferred with the same jump frequency data (see Fig. 9).In the case of N and C, only experimental bulk diffusion data are available in the literature.KMC data of N bulk diffusion are compared with those obtained using ion implantation and nuclear resonance in high-purity Ti samples (99.99%) 16 in Fig. 10.Additionally, C bulk diffusion evaluated with KMC is compared with experimental values obtained using diamond deposition on pure bulk Ti substrates 17 in Fig. 11.In both cases, the KMC simulation results are consistent with the experimental values.

A. Electronic structure of interstitial sites
It is a common belief that interstitial elements, such as O, N, and C, occupy octahedral and tetrahedral sites. 2 However, DFT calculations show that other interstitial sites can be TABLE V. Activation energies, Q (eV), of O, N, and C diffusion obtained using the MSD and KMC method.occupied by these elements.Our first-principles calculations confirm previous ab initio findings 11,32 that O occupies octahedral and less traditional crowdion and hexahedral sites.On the other hand, tetrahedral sites are unstable and relax to hexahedral sites.Based on the interstitial energies (see Table I), O shows a preference to occupy octahedral sites.The area under the total DOS curves in Fig. 12 in the anti-bonding region (above 0 eV) is much smaller compared to the pure Ti DOS when O occupies one of the interstitial sites, especially in the case of O occupying OC and HE sites.Additionally, the Fermi energy increases when O occupies OC and HE with respect to pure Ti (3.147 eV) as shown in Table VIII.
The DOS curve (Fig. 12) also reveals that the pseudo-gap at the Fermi energy is deeper when O sits at octahedral sites.These findings imply that Ti-O OC and Ti-O HE bonds are strong and stable with some covalent characteristics.On   16 sites, also preferring the octahedral sites amongst all possible interstitial sites (see Table I).In the anti-bonding region of the total DOS curves (Fig. 13), the areas under Ti-N OC and Ti-N HE curves are narrower than Ti-N BC and Ti-N CR .The pseudo-gap is deeper when nitrogen occupies hexahedral and octahedral sites than when it occupies crowdion and basal crowdion sites.The Fermi energy moves to the valence region in all cases, except for when nitrogen occupies the octahedral sites (Table VIII).Ti-N OC and Ti-N HE bonds therefore have a covalent nature, resulting in a stronger and more stable bond compared to when nitrogen occupies the BC and CR sites.Carbon also prefers to occupy octahedral sites like N and O; the other stable sites are crowdion and basal crowdion sites.The Ti-C OC bond is more covalent than the Ti-C CR and Ti-C BC bonds, as shown by total DOS curves in Fig. 14.The Ti-C OC DOS curve occupies a smaller portion of the anti-bonding region and its pseudo-gap is deeper.Nevertheless, C has some peculiarities when compared to N and O: the Fermi energies decrease compared to the pure Ti value (3.147 eV) when C occupies all possible interstitial sites, as shown in Table VIII.The highest energy level occupied by the Ti-C OC bond in the anti-bonding zone is larger than that of N and O. Also, the area under the Ti-C CR DOS curve in the anti-bonding region is smaller than the other two solutes.Therefore, the nature of Ti-C bond does not change between the different interstitial sites as much as it does for N and O. Ti-C OC bonds are weaker when compared to the Ti-O OC and Ti-N OC bonds.The detrimental influence of O, N, and C on ductility of Ti is well-documented and usually attributed to the crystal deformation induced by these solutes in the lattice. 2,33wasniak and co-workers 4 have suggested that the loss of plasticity is also a result of the electronic structure of solutesolvent bonds.They have studied the influence of O addition FIG. 13.Comparison of total DOS curves for pure Ti and supercells containing N at OC sites (labelled Ti-N OC ), HE sites (labelled Ti-N HE ), CR sites (labelled Ti-N CR ), and BC sites (labelled Ti-N BC ).The energies are referred to the Fermi energies (E F ) reported in Table VIII.VIII.
on the mechanical properties and ductility of α-Ti using ab initio, and their results show a strong relationship between Ti-O bond nature and plasticity.They find that a higher concentration of oxygen leads to a more covalent bond, reducing the ductility.The loss of due to additions of C is less dramatic when compared to N and O; this effect can be attributed to the more metallic nature of the Ti-C bonds. 33

B. Transition networks and anisotropy
The transition networks of O, N, and C diffusion in α-Ti were mapped using ab initio calculations and validated comparing KMC calculations of bulk diffusion with available experimental data.These three solutes show similar diffusion behaviours, significantly influenced by the stability of crowdion sites.The instability of tetrahedral sites significantly reduces the number of transition paths along c-axis.At the same time, the stability of crowdion sites allows jumps along The OC ↔ CR jump prevails over all the other jumps, as a result of its low energy barrier.In the case of C, the migration energy of CR → OC is significantly smaller than O and N, leading to markedly higher diffusivity values.On the other hand, the high energy barrier OC ↔ OC jumps along c-axis, due to the instability of basal octahedral sites, increases the activation energy for perpendicular diffusion.This leads to a very high diffusion anisotropy, as shown in Fig. 15, which is accentuated for C.

V. CONCLUSIONS
In conclusion, MSD and KMC methods informed by firstprinciples calculations are valid and valuable instruments to study and understand the interstitial migration mechanisms in α-Ti.The diffusivity values obtained with MSD and KMC are consistent with previous ab initio and experimental studies.First-principles calculations show that light elements tend to occupy octahedral sites and form semi-covalent bonds with Ti, which can lead to a reduction of plasticity in Ti alloys.Also, less traditional interstitial sites are stable over tetrahedral sites: hexahedral, crowdion, and basal crowdion sites can all be occupied by these solute elements.Possible pathway networks have been mapped between these interstitial sites, showing a similar diffusion behaviour of O, N, and C. The diffusion of these elements is mainly controlled by the jump occurring between OC and CR sites within the basal plane, resulting in high diffusion anisotropy.The transition networks of these solutes lay the foundations for the control and prevention of α-case formation and growth in Ti alloys.
FIG. 3. Schematic representation of stable interstitial sites occupied by (a) O, (b) N, and (c) C in α-Ti unit cell.

FIG. 14 .
FIG.14.Comparison of total DOS curves for pure Ti and supercells containing C at OC sites (labelled Ti-C OC ), CR sites (labelled Ti-C CR ), and BC sites (labelled Ti-C BC ).The energies are referred to the Fermi energies (E F ) reported in TableVIII.

TABLE I .
Interstitial energies, E i (eV), of O, N, and C in α-Ti.The notation in brackets indicates the interstitial position occupied after the relaxation.
) and 4(d)-4(f).When C occupies OC, it can perform two OC → OC, six OC → CR, or six OC → BC jumps.From CR sites, two CR → OC and four

TABLE III .
Migration energy, ∆H m (eV), and attempt frequency, ν (THz), of N interstitial jumps in α-Ti.BC jumps are possible.Finally, from BC sites, four BC → CR and two BC → OC jumps are possible.TableIVshows the migration energies and attempt frequencies of these jumps.

TABLE IV .
15gration energy, ∆H m (eV), and attempt frequency, ν (THz), of C interstitial jumps in α-Ti.haveexperimentallyevaluated the bulk diffusivity using ion implantation and nuclear resonance in high-purity Ti samples (99.99%) annealed for 10 days at 1133 K in order to be able to neglect the grain-boundary diffusion and residual stress effect.15Theyhaveextrapolated the following bulk diffusion parameters: Q = 1.752 eV and D 0 = 2 × 10 −7 m 2 /s.As shown in Figs.5 and 8, our findings are consistent with these two studies.In particular, the KMC method predicts the bulk diffusion very well.Our data match the results of Wu and Trinkle in the case of perpendicular diffusion.On the other hand, the parallel diffusivity values evaluated in this work overestimate Wu and Trinkle's values.The migration energies for CR → OC and CR → HE jumps are much smaller in the study of Wu and Trinkle: 0.28 eV and 0.24 eV, respectively.

TABLE VII .
Comparison of migration energies of O interstitial transitions in α-Ti computed with (i) 54-atom supercells, 350 eV cutoff energy, and constant-shape relaxation (∆H cs m ), (ii) 96-atom supercells, 400 eV cutoff energy, and fully relaxed supercells, (∆H), with the results of Wu and Trinkle, 11 ∆H Wu-USPP . 9. Parel diffusivity of O evaluated using Wu and Trinkle's formulation (D Formula ) and KMC simulations (D KMC ) neglecting CR → OC and CR → HE jumps both inferred with the jump frequency values evaluated in this work, and with Wu and Trinkle's jump frequencies data (D Formula FIG FIG. 12.Comparison of total DOS curves for pure Ti and supercells containing O at OC sites (labelled Ti-O OC ), HE sites (labelled Ti-O HE ), and CR sites (labelled Ti-O CR ).The energies are referred to the Fermi energies (E F ) reported in Table VIII.TABLE VIII.Fermi energy (eV) of supercells containing O, N and C at stable interstitial sites in α-Ti.