Light element diffusion in Mg using first principles calculations: Anisotropy and elastodiffusion
LLight element diffusion in Mg using first principles calculations:Anisotropy and elastodiffusion
Ravi Agarwal and Dallas R. Trinkle ∗ Department of Materials Science and Engineering,University of Illinois at Urbana-Champaign, Urbana, IL 61801 (Dated: September 20, 2018)
Abstract
The light elemental solutes B, C, N, and O can penetrate the surface of Mg alloys and dif-fuse during heat treatment or high temperature application, forming undesirable compounds. Weinvestigate the diffusion of these solutes by determining their stable interstitial sites and the inter-penetrating network formed by these sites. We use density functional theory (DFT) to calculatethe site energies, migration barriers, and attempt frequencies for these networks to inform ouranalytical model for bulk diffusion. Due to the nature of the networks, O diffuses isotropically,while B, C, and N diffuse anisotropically. We compute the elastodiffusion tensor which quantifieschanges in diffusivity due to small strains that perturb the diffusion network geometry and themigration barriers. The DFT-computed elastic dipole tensor which quantifies the change in siteenergies and migration barriers due to small strains is used as an input to determine the elasto-diffusion tensor. We employ the elastodiffusion tensor to determine the effect of thermal strains oninterstitial diffusion and find that B, C, and N diffusivity increases on crystal expansion, while Odiffusivity decreases. From the elastodiffusion and compliance tensors we calculate the activationvolume of diffusion and find that it is positive and anisotropic for B, C and N diffusion, whereas itis negative and isotropic for O diffusion.
PACS numbers: 66.10.C-, 66.10.cg, 66.30.J-, 66.30.Ny a r X i v : . [ c ond - m a t . m t r l - s c i ] J u l . INTRODUCTION Magnesium and its alloys have found increased application in the automotive industry dueto their higher strength-to-weight ratio than steel and aluminum alloys, which reduces vehicleweight leading to increase in fuel efficiency . Mg alloys interact with the surroundinggaseous atmosphere during their application which can lead to the penetration of lightimpurity elements. These impurities can also get introduced due to interaction with reactivegases during heat treatment, leading to the formation of oxide layers on the surface orprecipitates at grain boundaries which can be detrimental to strength . Experiments haveshown that O, C and N can react with Mg to form oxides, carbides and nitrides . Boron isused for Fe removal during Mg processing , but a small amount of B may be retained as animpurity. The penetration of these impurities into bulk is governed by thermally activatedprocesses and a detailed study of their diffusion mechanisms can provide insights that mayhelp to mitigate them.There have been few theoretical studies on the behavior of light elements in hcp metals.Wu et al. studied the influence of substitutional B, C, N and O on the stacking faults andsurfaces of Mg using density functional theory (DFT). All four elements reduce the unstablestacking fault energy and surface energy of Mg and enhance the ductility according to theRice criterion, with O having the largest impact . Atomisitic studies of light elements inhcp metals—O in α -Ti , O and N diffusion in α -Hf , and O in multiple hcp metals —modeled the diffusion of solutes through the networks formed by interstitial sites. However,a theoretical or experimental study of interstitial diffusion in Mg is absent except for thelimited experimental data for C diffusion .We analyze the diffusion of B, C, N and O in the dilute limit in hcp Mg using DFTcalculations to inform an analytical diffusion model . We also study the changes indiffusivities due to strain from thermal expansion. Section II details the DFT parametersused to determine the energetics of interstitial sites and the migration barriers betweenthem. Section III lays out the inputs for the diffusion model: probabilities of occupyingsites, connectivity networks between these sites and the transition rates for these networks.We derive analytical expressions for interstitial diffusivity in hcp crystals and apply them todiffusion of B, C, N and O in Mg. We find that the O diffusion is isotropic while B, C, and Ndiffusion is anisotropic. Section IV discusses the elastic dipole tensors of solutes at interstitial2ites and transition states, which determine the changes in the transition energetics of solutesdue to small strains. Section V defines the elastodiffusion tensor , which quantifies theeffect of small strains on diffusivity and discusses the sign inversion behavior of elastodiffusioncomponents with temperature. We find that the activation volume of O diffusion is negativewhich leads to an increase in O diffusion under hydrostatic pressure. We also find thatthe diffusivity of O decreases with thermal expansion while the diffusivity of B, C and Nincreases. II. COMPUTATIONAL DETAILS
We perform the DFT calculations using the Vienna ab-initio simulation package vasp which is based on plane wave basis sets. The projector-augmented wave psuedopotentials generated by Kresse describe the nuclei and the valence electrons of solutes and Mg atoms.The solute atoms B, C, N, and O are described by [He] core with 3, 4, 5 and 6 valenceelectrons respectively. We use the [Ne] core with 2 valence electrons for Mg instead ofthe [Be] core with 8 valence electrons because the energies computed using either choiceof psuedopotential differ by less than 20 meV. Electron exchange and correlation is treatedusing the PBE generalized gradient approximation. We use a 4 × × × × k -point mesh to sample the Brillouin zone.Methfessel-Paxton smearing is used with energy width of 0.25 eV to integrate the densityof states; the k-point density and smearing width are based on convergence of the DOScompared with tetrahedron integration. A plane wave energy cutoff of 500 eV is requiredto give an energy convergence of less than 1 meV/atom. All the atoms are relaxed using aconjugate gradient method until each force is less than 5 meV/˚A. The Mg unit cell has ahexagonal close packed (HCP) crystal structure with DFT calculated lattice parameters of a = 3 .
189 ˚A and c/a ratio of 1.627 which agree well with values reported from experiments, a = 3 .
19 ˚A and c/a = 1 . .We use DFT to calculate the energy of solutes at various sites and use the climbing-imagenudged elastic band (CNEB) method to locate the transition states between the sites. Thesite (or solution) energy E α of a solute X at an interstitial site α is the difference betweenthe energy of a Mg supercell containing solute X at site α , E (Mg + X α ), and the energy3f a pure Mg supercell, E (Mg ), E α = E (Mg + X α ) − E (Mg ) . (1)We also determine the site energy for a solute X as a substitutional defect, E sub , E sub = E (Mg + X sub1 ) − E (Mg ) (2)where E (Mg + X sub1 ) is the energy of supercell where one of the Mg atoms is substituted bya solute atom X. Both the interstitial site energy E α and the substitutional site energy E sub for solute X are referenced to its elemental state. The energy differences ∆ E = E α − E sub for the solutes B, C, N and O are –1.48, –3.23, –4.34 and –4.19 eV, where α is the interstitialsite with the lowest energy, and is independent of the reference state for the solutes. Since,the energies of interstitial sites are lower than the substitutional site, these solutes are likelyto diffuse through networks of interstitial sites. We use CNEB with one image to locatethe transition state between two interstitial sites. Similar to Eq. 1, the energy E α - β of thetransition state between site α to site β is referenced to the elemental state of X E α - β = E α - β (Mg + X ) − E (Mg ) (3)where E α - β (Mg + X ) is the energy at the transition state obtained from a CNEB calcula-tion. We report the interstitial site energies and the transition state energies relative to theinterstitial site with the lowest energy, which is independent of the reference state for thesolutes. III. DIFFUSION MODEL
We calculate the occupation probabilities at interstitial sites and transition rates for dif-fusion pathways from DFT-computed site energies, transition state energies, and vibrationalfrequencies. The probability ρ α of a solute occupying a particular site α at temperature T is ρ α = ν ∗ α · exp( − E α /k B T ) (cid:80) β ν ∗ β · exp( − E β /k B T ) , (4)where k B is the Boltzmann constant, in the denominator is the normalization constantsummed over all the interstitial sites in the unit cell and ν ∗ α is the site prefactor proportional4o the Arrhenius factor for formation entropy of site α , exp ( S α /k B ), calculated from thevibrational frequencies ν ∗ α = 1 (cid:81) p =1 ν α,p . (5)This expression ignores interstitial-interstitial interaction, and is exact in the dilute con-centration limit. We compute the vibrational frequencies of a state using the one atomapproximation by diagonalizing the dynamical matrices corresponding to the interstitialatom. The dynamical matrices are obtained from the forces induced on interstitial atomsby small displacements ( ± .
01 ˚A) from their equilibrium positions, while keeping the otheratoms fixed. From transition state theory, the rate λ α - β for a solute to transition from site α to site β at temperature T is λ α - β = ν ∗ α - β · exp( − ( E α - β − E α ) /k B T ) . (6)The attempt frequency ν ∗ α - β for the α to β transition is calculated using the Vineyardexpression , which is the product of vibrational frequencies ν α,p at the initial site α di-vided by the product of real vibrational frequencies ν α - β,q at the transition state ν ∗ α - β = (cid:81) p =1 ν α,p (cid:81) q =1 ν α - β,q . (7)At equilibrium, the transition between site α and site β obeys detailed balance ρ α · λ α - β = ρ β · λ β - α . (8)Figure 1 shows the newly found distorted hexahedral dh site in Mg along with the otherinterstitial sites (h, t, c, o) which have been discussed previously for O in α -Ti . The dhsite is stable for B and C, and is located between two nearest Mg atoms in the basal planewith a displacement of 0.17 ˚A for B and 0.40 ˚A for C towards the nearest hexahedral hsite. The h site has three basal Mg neighbors and two other Mg neighbors located directlyabove and below it, which are further away. The four-atom coordinated tetrahedral t siteis stable for O and lies 0.65 ˚A along the c direction from an basal plane containing three ofits Mg neighbors. The six-atom coordinated octahedral o site is stable for all four solutes.The six-atom coordinated non-basal crowdion c with lower symmetry than o site has twonearest neighboring Mg atoms lying in adjacent basal planes which get displaced away fromthe c site while the other four neighbors lying further apart get displaced towards the c siteon relaxation. The c site is stable for C and N but unstable for B and O.5 IG. 1. (color online) Positions of interstitial sites in the unit cell of hcp Mg. The octahedral (o,orange), tetrahedral (t, red), hexahedral (h, blue), distorted hexahedral (dh, cyan), and crowdion(c, yellow) interstitial sites are shown relative to host Mg atoms (Mg, white). In an hcp unit cell,there are two o, two h, four t, six c and six dh sites. The transitions between stable interstitialsites determine the possible diffusion pathways. The unit cell vectors a and a form the basalplane (0001) and the vector c (also referred as the c -axis ) is perpendicular to it. Figure 2 shows the possible diffusion networks between interstitial sites for hcp systems,which are inputs to our diffusion model . A solute at a o site can jump to the followingneighboring sites: two o sites lying above and below along the c -axis with transition rate λ o-o( c ) ; six o sites lying in the same basal plane with λ o-o(b) in cases where the c site isunstable ; six neighboring c sites with λ o-c ; six h sites with λ o-h ; six t sites with λ o-t and sixdh sites with λ o-dh . A solute at a h site can jump to: six o sites with λ h-o ; six c sites with λ h-c and three dh sites lying in the same basal plane with λ h-dh . The c site is between twoh sites which lie in adjacent basal planes and also between two o sites in the same basalplane. A solute from a c site can jump to those neighboring o and h sites with λ c-o and λ c-h .A solute at a t site can jump to three neighboring o sites which are all lying either aboveor below the t site with λ t-o , and to one neighboring t site lying either above or below with λ t-t . A solute at a dh site can jump to one neighboring h site with λ dh-h and to two nearestdh sites in the same basal plane with λ dh-dh .Figure 3 shows the energies for the interstitial sites and the transition states of activediffusion pathways for all four solutes. Active diffusion pathways for a solute are determinedby its set of stable sites. The set of stable sites for B is { o, dh } , for C it is { o, h, c, dh } , forN it is { o, h, c } and for O it is { o, t } . All DFT energies are relative to the lowest-energysite which is the ground state. The o site is the ground state for B, C and N, while the t6 -o( c ) o-o(b), o-c h-c o-ht-t, t-o o-dh, dh-h o-dh, dh-dhFIG. 2. (color online) Interstitial sites and site-to-site connectivity in hcp crystals. Connections be-tween two neighboring sites form diffusion pathways which are shown as lines colored correspondingto the colors of the interstitial sites. For example, o-o( c ) and o-o(b) are octahedral site-to-octahedralsite diffusion pathways along the c -axis and in the basal plane of hcp Mg, respectively. The dif-fusion pathways shown in the top row are un-correlated, while correlated diffusion pathways areshown in the bottom row. These correlated pathways are the combined connections formed amongo and t sites, o, dh and h sites, and o and dh sites. In the last two figures of the bottom row, the c -axis is tilted and the cell is rotated counter-clockwise around the c -axis for better visibility ofconnections and sites. site is the ground state for O. The transition between two sites is shown as a line connectionand the associated value is the transition state energy. For example, in the case of O, t isthe ground state and o is metastable with energy 0.21 eV. The active diffusion pathways forO (refer to Fig. 2) are o-o, t-t (both along the c -axis), and t-o with transition state energiesof 1.01, 0.09 and 0.7 eV respectively. Since there is no direct o-o (b) jump in the basalplane—which would pass through the unstable c site—basal diffusion occurs by combiningo-t and t-o jumps.Table I lists analytical expressions for diffusivity based on the active diffusion pathways7 .00.51.01.5 Boron o dh c ) E (eV) Carbon o h c dh1.65 E (eV) Nitrogen o h c1.64 E (eV) Oxygen o t0.09 E (eV) FIG. 3. Energetics of stable sites and the transition states between them, relative to the lowest-energy interstitial site for B, C, N, and O solutes in Mg. Interstitial sites are marked on thehorizontal axis, and their relative site energies are shown in bold below the thick horizontal baselines. Thin lines from one site to another (or the same) site denote transitions, and the associatednumber is the energy at the transition state between those two sites. For example, in the case ofB, the o site is the lowest energy site and the energy of the metastable dh site relative to it is 0.90eV. Thin lines starting and ending from the thick base line of o denotes the o-o transition. Theassociated transition state energies in eV are 1.08( c ) for the transition along the c -axis and 0.73(b)for the transition in the basal plane. formed by the stable sites, in terms of occupation probabilities and transition rates. Wefollow the approach of near-equilibrium thermodynamics to calculate the diffusivity D byfinding a steady state solution for the system in equilibrium distribution with a small pertur-bation in the chemical potential gradient of the solute . The derived analytical expressionsfor solute diffusivity are made up of bare mobilities and correlation effects. Table I lists theterm-by-term contributions to the basal diffusivity D b and the c -axis diffusivity D c from eachtype of transition. The bare mobility terms have the form of a site probability multipliedby a transition rate. The correlation effects are present in dh-o, dh-dh and dh-h transitionswhich contribute to the basal diffusivity as well as in t-o and t-t transitions which con-tribute to the c -axis diffusivity. Each of these networks show correlation as the jumps from8 ABLE I. Analytical expressions for interstitial solute diffusivity in the basal plane ( D b ) and alongthe c -axis ( D c ) through the network formed by interstitial sites in the hcp crystal. These expressionsare functions of transition rates ( λ ) between interstitial sites and the occupation probability of eachtype of interstitial site. The occupation probability of each type of site is the product of ρ (fromEq. 4) and its multiplicity in the unit cell. The occupation probability for any o, h, t, dh and c siteis 2 ρ o , 2 ρ h , 4 ρ t , 6 ρ dh and 6 ρ c , respectively. These analytical expressions for diffusivity are validfor any interstitial solute diffusing in an hcp crystal with lattice parameters a and c and having aset of stable interstitial sites corresponding with that network for a Markovian diffusion process.network a − · D b c − · D c o, dh 2 ρ o λ o-o(b) ρ o λ o-dh λ dh-dh λ dh-o + 3 λ dh-dh ρ o λ o-o( c ) ρ o λ o-dh ρ o λ o-c ρ h λ h-c ρ o λ o-dh λ dh-h λ dh-o + λ dh-h ρ o λ o-o( c ) ρ o λ o-dh ρ h λ h-c ρ o λ o-c ρ h λ h-c ρ o λ o-h ρ o λ o-o( c ) ρ o λ o-h ρ h λ h-c ρ t λ t-o ρ o λ o-o( c ) ρ t λ t-o λ t-t λ t-o + 16 λ t-t particular sites (dh and t) are unbalanced : the sum (cid:80) β λ α - β δ x α − β (cid:54) = 0 for displacements δ x α − β from site α to β . This leads to a correlated random walk where, for example, if aninterstitial is in a tetrahedral site with a low t-t barrier it is very likely to be in that sametetrahedral site after two jumps; hence, a large (anti)correlation between the displacementvector in subsequent jumps. The analytical expressions are applicable in any hcp crystal forany solute having a set of interstitial sites corresponding with that network for a Markoviandiffusion process. Our expression for the set of sites { o, h, c } agree with the expression for Odiffusing in α -Ti . In the case of t-t jumps which tend to have low barriers, the assumptionof “independent” tetrahedral sites becomes invalid; instead, the pair is similar to a super-basin which thermalizes rapidly, and the λ t-t disappears from the diffusivity as λ t-t → ∞ .The site energies and site prefactors, as well as the attempt frequencies and transition stateenergies of all the transitions for B, C, N and O, is available in tabular form .Figure 4 shows that O diffuses isotropically while B, C, and N diffuse anisotropically. Band C diffuse faster in the basal plane than along the c -axis while N diffuses faster alongthe c -axis than in the basal plane. The analytical expressions in Table I give the diffusivity9 /T . − − − − − − − D P V − D Bb D Cb D Nb D Ob D B c D C c D N c D O c T . D CExperiment
FIG. 4. (color online) Analytical results for diffusivities in the basal plane ( D Xb ) and along the c -axis( D X c ) of Mg for interstitial solute X = B, C, N and O. Diffusion of O is isotropic while diffusion of Band C is slower along the c -axis than in the basal plane and diffusion of N is faster along the c -axisthan in the basal plane. The analytical expressions listed in Table I are employed to computethe variation of diffusivity with temperature. Also shown is the diffusivity of C ( D CExperiment ),determined experimentally by Zotov et.al at four temperatures between 773–873K. as a function of temperature. For all temperatures from 300K to 923K (the melting pointof Mg), the basal diffusivities of the four solutes follow D Bb > D Ob > D Nb ≈ D Cb and the c -axis diffusivities follow D O c > D B c > D N c > D C c . Zotov et.al measured the diffusivityof C experimentally in the temperature range of 773–873K (500–600 ◦ C) and our resultsoverestimate their measured diffusivity by a factor of 10–80. With only the single experimentfor comparison, it is difficult to assess the source of the discrepancy.Table II lists the activation energies and diffusivity prefactors obtained from Arrheniusfits to the diffusivity plots (Fig. 4). For each solute, the comparison between the activationenergy for diffusion Q and the migration energies of individual transitions (see Fig. 3) indi-cates the dominant type of transition that contributes most to diffusion. In the case of O,the migration energy of t-o transition is 0.70 eV which is close to the activation energy of0.69 eV, so this transition contributes more than the other transitions to both diffusivities.Similarly, o-o basal and o-c transitions dominate for basal diffusion of B and C, respectively,while o-dh transitions dominate for c -axis diffusion of both these solutes. However, for N,all transitions except o-o along c axis have similar energies, so it is likely that more than10 ABLE II. The Arrhenius fitting parameters for basal ( D Xb ) and c -axis ( D X c ) diffusivities throughactive networks of sites for interstitial solute X = B, C, N, and O. The diffusivities vary withtemperature according to the Arrhenius model D = D · exp( − Q/k B T ), where D is the diffusivityprefactor, Q is the activation energy of diffusion, T is temperature in K, and k B is the Boltzmannconstant. The comparison of energy barriers from Fig. 3 to the activation energy Q gives thedominant transition.Solute Network D Xb D X c X D (m s − ) Q (eV) D (m s − ) Q (eV)B o, dh 2.52 × − × − × − × − × − × − × − × − one transition type contributes to both diffusivities. IV. ELASTIC DIPOLE TENSOR
The elastic dipole tensor quantifies the elastic interaction energy between an externalstrain field and the point defect in the small strain limit. The dipole tensor is equal tothe negative derivative of elastic energy E with respect to strain ε . The elastic dipolecomponents P ij are computed from the stress tensor σ after relaxing the ions while keepingthe supercell shape and volume V fixed in the presence of the interstitial , P ij = − dEdε ij ≈ σ ij V. (9)The elastic dipole tensor determines the change in site energies and transition state ener-gies of interstitial solutes due to small strain. The site energy E α ( s ) ( ε ) of α with orientationvector s under small strain ε is approximated by the linear relation E α ( s ) ( ε ) ≈ E α (0) − (cid:88) ij P α ( s ) ,ij ε ij , (10)where E α (0) is the site energy of α in the unstrained cell and P α ( s ) ,ij are the elastic dipolecomponents of site α with orientation s . In the infinitesimal strain limit, the sites and11etwork topology remains unchanged; with larger finite strains, sites may become unstableor change the network topology, which requires a new analysis of network. The vector s distinguishes the multiple sites of the same type which are present in an hcp unit cell. Theorientation of c site is defined as the vector connecting it to the nearest o site and theorientation of dh site is defined as the vector connecting it to the nearest h site. In a hcpunit cell (see Fig. 1), there are two o, two h, four t, six c and six dh sites. In an unstrainedcell, multiple sites of the same type have the same energy. However, strain can cause thesesites to become nonequivalent in energy depending on their elastic dipole tensor which maydepend on their site orientation. The dipoles for o, h, and t sites are independent of theirorientation vector while the dipole for c and dh sites depend on their orientation vector.Similarly, the transition state energy E v α ( s )- β ( s (cid:48) ) ( ε ) for site α of orientation s to site β oforientation s (cid:48) under strain is E v α ( s )- β ( s (cid:48) ) ( ε ) ≈ E α - β (0) − (cid:88) ij P v α ( s )- β ( s (cid:48) ) ,ij ε ij , (11)where v is the vector from site α to β , E α - β (0) is the v -independent transition state energyin the unstrained cell and P v α ( s )- β ( s (cid:48) ) ,ij are the elastic dipole components at the transitionstate corresponding to v . As discussed previously in Fig. 2, there are multiple transitionsof the same type distinguished through their transition vectors v . In a strained cell, thesetransitions can have different transition state energies depending on their dipole tensorswhich may depend on their transition vectors.Tables III and IV list the components of the elastic dipole tensor at representative inter-stitial sites with orientations s , and representative transition states with transition vectors v .We diagonalize the elastic dipole tensors along three principal axes ( e , e , e ), and reportthe diagonalized entries entries ( P , P , P ) and principal axes. From Table III, the elasticdipole components in the two orthogonal basal directions are equal for o, h, and t sites dueto the basal symmetry of these sites. The trace of the elastic dipole for N and O at o sites isnegative, leading to the volumetric contraction upon cell relaxation, in contrast to the otherinterstitial sites. The ground state configuration of N undergoes volume contraction on cellrelaxation while the ground state configuration of B, C, and O undergoes volume expansionon cell relaxation. In the case of the dh site, its two nearest Mg atoms experience largeratomic forces compared to other Mg atoms therefore, the elastic dipole for the dh site hasthe largest component in the [1120] direction which connects these two nearest Mg atoms.12 ABLE III. Elastic dipole tensors P at representative interstitial sites for B, C, N, and O in Mg.The symmetric elastic dipole tensor is diagonal along three principal axes e , e , and e and hasunits of eV. For c and dh sites, the dipole tensors and their axes depend on the orientations s of thesites with respect to the nearest o and h sites, respectively, whereas the dipole tensors for o, t andh sites are independent of orientation. The possible orientations of dh sites with respect to an hsite are [1100], [1010] and [0110], and the orientations of c sites with respect to an o site are [2110],[1120] and [1210]. Here the dipole tensor of each type of site is given for one representative s , andother tensors with different s are obtained by applying the appropriate point group operations onthe representative dipole tensor.Solute Site Orientation ( s ) P P P e e e B o any 2.38 2.38 2.55 orthogonal basal vectors [0001]dh [1100] 11.03 0.04 − .
49 [1120] [1100] [0001]C o any 1.08 1.08 0.24 orthogonal basal vectors [0001]h any 4.74 4.74 − .
10 orthogonal basal vectors [0001]c [2110] 6.59 4.20 − .
18 [0
13 13 12 ] [2110] [011 ]dh [1100] 8.94 − . − .
86 [1120] [1100] [0001]N o any 0.00 0.00 − .
39 orthogonal basal vectors [0001]h any 3.22 3.22 − .
81 orthogonal basal vectors [0001]c [2110] 4.22 4.31 − .
39 [0
13 13 12 ] [2110] [011 ]O o any − . − . − .
76 orthogonal basal vectors [0001]t any 2.06 2.06 0.79 orthogonal basal vectors [0001]
From Table IV, most of the transition states break the symmetry of the crystal except forthe o-o and t-t transitions along the c -axis which obey the basal symmetry. Because ofthe basal symmetry, the transition state energies of the o-o ( c -axis) and t-t transitions withdifferent v remain equivalent in the strained cell while the same is not true for the othertypes of transitions.Elastic dipole tensors for symmetry-equivalent sites with different s , and symmetry-equivalent transitions with different v , are obtained by point group operations on the rep-resentative dipole tensors in Tables III and IV. For example, the three c sites in the basalplane with different orientations ([2110], [1120] and [1210]) are all related Wyckoff sites, that13 ABLE IV. Elastic dipole tensors P at representative transition states for B, C, N, and O in Mg.The transition state from site α to site β is denoted by α - β , and v is the vector connecting thesetwo sites. The symmetric elastic dipole tensor is diagonal along three principal axes e , e , and e and has units of eV. The dipole tensor of an equivalent transition with a different v is obtainedby applying the appropriate point group operation to the given dipole tensor. The variable x forB and C is 0.197 and 0.238, respectively, and variable z for O is 0.153. The values of x and z areobtained from the relaxed position of dh and t sites in the Mg supercell, respectively.Solute α - β Transition ( v ) P P P e e e B o-o [000 ] 5 .
34 5 . − .
58 orthogonal basal vectors [0001]o-o [2110] − .
08 2 .
74 8 .
56 [011 ] [2110] [0
13 13 12 ]o-dh [ xx ] 10.99 0.19 − .
61 [1120] [110 ] [110 ]dh-dh ( x − )[2110] 7 .
69 4 . − .
20 [0110] [2110] [0001]C o-o [000 ] 3 .
63 3 . − .
22 orthogonal basal vectors [0001]o-c [2110]] − .
93 1 .
57 5 .
19 [011 ] [2110] [0
13 13 12 ]o-dh [ xx ] − .
60 0 .
58 7 .
28 [110 (cid:113) ] [110 (cid:113) ] [1120]h-c [0
16 16 14 ] 3 . − .
68 4 .
45 [2110] [0
16 16 14 ] [011 ]h-dh ( x − )[1100] 7 .
59 1 . − .
01 [1120] [1100] [0001]N o-o [000 ] 2 .
58 2 . − .
16 orthogonal basal vectors [0001]o-c [2110]] − .
08 3 .
23 3 .
58 [011 ] [2110] [0
13 13 12 ]o-h [
13 13 ] − .
08 0 .
31 4 .
61 [0 . , . , , .
10] [0 . , . , , .
60] [1120]h-c [0
16 16 14 ] 3 . − .
14 3 .
94 [2110] [011 ] [0
13 13 12 ]O o-o [000 ] 2 .
37 2 .
37 1 .
76 orthogonal basal vectors [0001]t-t [000( − z )] 1 .
97 1 . − .
67 orthogonal basal vectors [0001]o-t [
13 13 z ] 0 .
67 1 .
37 2 .
07 [0 . , . , , .
60] [0 . , . , , .
12] [1120] are transformed by 120 ◦ rotations about the c -axis; call that transformation matrix R . Thedipole tensors for the other two equivalent sites s (cid:48) are P α ( s (cid:48) ) = RP α ( s ) R T (12)where P α ( s ) is the representative dipole tensor and R transforms s to s (cid:48) . Similarly, the dipole14ensors of all the other sites are calculated using their associated transformation matrices.The same operations are carried out for all the transition state dipole tensors based onthe symmetry of the transition vectors v . The dipole data in Cartesian basis for all theseequivalent sites and equivalent transitions for B, C, N and O are available in tabular form .This dipole tensor data is used to estimate changes in site energies and the changes inmigration barriers of transitions under strain using Eqs. 10 and 11, which are inputs to theelastodiffusion tensor calculations. V. ELASTODIFFUSION TENSOR
Strain affects the diffusivity of solutes by changing the jump vectors and migration bar-riers of the diffusion network. The first order strain dependence of diffusivity is representedwith the elastodiffusion tensor d d ijkl = ∂D ij ∂ε kl , (13)and is derived using perturbation theory . The contribution d geom to the elastodiffusiontensor from the changes in jump vectors is d geom ijkl = 12 ( D jk (0) δ il + D il (0) δ jk + D ik (0) δ jl + D jl (0) δ ik ) , (14)where δ ij are the Kronecker deltas. Hence, if the diffusivity has Arrhenius temperaturedependence, then so does the geometric term in the elastodiffusion tensor. The contribution d mb from changes in the migration barriers is determined by the elastic dipole tensors ofthe migration barriers and sites. The elastic dipole tensor of a transition state relative toinitial site determines the rate of that transition under strain and the elastic dipole tensorof interstitial site determine the occupation probability of that site under strain. The term d mb is the sum of contributions from each transition; these contributions are proportionalto the product of the inverse temperature, transition rate, and difference of transition statedipole and thermal average dipole of interstitial sites. The contribution from one transitioncan be represented as d k B T · exp( − E/k B T ) (15)where the elastic dipole terms are absorbed in the “prefactor” d , which has units of eV · m s − , and E is the barrier of the dominant transition.15 ABLE V. The fitting parameters d and E in Eqn. 15 for the components of the B, C, N, andO elastodiffusion tensor in Mg over 300–923K. The elastodiffusion tensor in Voigt notation has sixunique components in an hcp crystal, where d = ( d − d ) /
2. A subset of components changesign with temperature; their transition temperature is listed in lieu of fitting parameters (c.f., Fig. 5for the temperature dependence). The “activation barrier” E corresponds closely to the migrationbarrier of the dominant transition. The d and d components for B, all diagonal components forC, and d and d for N and O are negative throughout the temperature range (i.e. have negative d ). The negative d implies that the increase in diffusivity caused by lowered migration barriersis greater than the decrease in diffusivity due to reduced jump vectors under compressive strains.For d of B, the geometric contribution is dominant and is best described with an Arrhenius fitof 1 . × − m / s · exp( − . /k B T ).B C N O d (eVm s − ) E (eV) d (eVm s − ) E (eV) d (eVm s − ) E (eV) d (eVm s − ) E (eV) d (854.7K) − . × − d − . × − . × − . × − d . × − . × − . × − d . × − . × − . × − d − . × − − . × − . × − d . × − ∗ − . × − − . × − − . × − d . × − − . × − − . × − − . × − The symmetry of the hexagonal closed-packed crystal reduce the number of unique elasto-diffusion components to six. We use Voigt notation, similar to elastic constants, to representthe indices of the fourth rank tensor as both diffusivity and strain are symmetric secondrank tensors. The reduction by symmetry is the same as the elastic constants, except that d ij is not necessarily equal to d ji . In the case of hcp, the non-zero elastodiffusion elementsare d = d , d , d , d = d , d = d , d = d , and d = ( d − d ) /
2. The changein jump vectors contributes only to d , d , d , and d . Unlike the contribution from thechange in jump vectors, the change in migration barrier can contribute to all six independentcomponents of elastodiffusion tensor and need not only be positive.Table V shows that the contribution d mb dominates over the contribution d geom due to the16elatively larger values of elastic dipole tensor components compared to k B T (see eqn. 15) forall the temperatures between 300–923K. However, the contribution d geom is greater than thecontribution d mb for the d component for B due to larger transition rate of o-o transitionin basal plane and for the d component for B and O at temperatures above crossover(discussed in the later paragraph). Equation 15 is used to fit the elastodiffusion componentbecause of the larger contribution from d mb over d geom and also due to the dominant transitionfor each solute. The fitting parameter E in Table V corresponds to the migration barrierof the dominant transition. These dominant transitions under strain is same as that in theunstrained crystal, except for the basal components d , d and d for C which are nowdominated by the h-dh transition. The remaining basal component d of C is governed byo-c transition and the basal components ( d , d and d ) and d of B are governed byo-o(b)transition. The non-basal components ( d and d ) are governed by o-dh transitionfor both B and C. The isotropic o-t transition is dominant for all the components for O,and in N, both o-h and h-c transitions, which have similar migration barriers, contribute toelastodiffusion components.Figure 5 shows that five of the elastodiffusion components for oxygen change sign (fewerfor B, and N) due to the small energy separation from the ground state and the metastablestates, while for B, C, and N the energy separation is significant. The change in sign frompositive (filled symbol) to negative (unfilled symbol) is observed as dips in the logarithmof the magnitude of d and the associated crossover temperature is listed in parenthesisin Table V for these components. The sign inversion of these components is due to thecompeting mechanism dominating over different temperature which we observe as differentslopes on opposite side of crossover. The sign inversion of d , d , d , and d for O is due tothe large variation in thermally averaged elastic dipole tensor of sites, which occurs becauseof the low energy separation of 0.21 eV between o and t sites. The difference between thetransition state dipole and the thermally averaged dipole contributes to the elastodiffusioncomponent sign changes with temperature as the o and t sites have different elastic dipoles.However, for d for B and O, the sign inversion is due to the competition between thenegative contribution of d mb and positive contribution of d geom , where the former dominatesbelow the crossover temperature (due to smaller value of k B T compare to dipole tensor, c.f.Eqn. 15) and the latter dominates above the crossover temperature. For the component d of N, sign inversion is due to the o-c transition dominating above the crossover temperature17 d d d d positive negative FIG. 5. (color online) Components of the elastodiffusion tensor d that change sign as a function oftemperature, for B, N, and O. The magnitudes of each component are shown with filled symbolsfor positive values and unfilled for negative values. For a component, changes of sign is observedas a dip in the curve and the crossover temperatures is listed in Table V. The sign inversionof these components is caused by two competing mechanisms, which dominate at either low orhigh temperatures. Five components of the elastodiffusion tensor for O change sign, and eachcomponent has a different crossover temperature. while the o-h transition dominates below the crossover. The sign inversion behavior ofdifferent components suggest that the diffusivity under strain will have contrasting featuresaround a specific temperature which we observe for the activation volume of diffusion andfor the effect of thermal expansion on diffusion. Activation volume of diffusion
The elastodiffusion tensor together with the elastic compliance tensor computes the ac-tivation volume of diffusion. The activation volume of diffusion V ij describes the pressure p dependence of diffusivity as D ij ( p ) = D ij (0) · exp( − pV ij k B T ) , (16)18 /T . $ F W L Y D W L R Q Y R O X P H Ω V Bb V B c V Cb V C c V Nb V N c V Ob V O c T . FIG. 6. (color online) Activation volume for basal diffusion V b and c -axis diffusion V c , relative tothe Mg atomic volume Ω = 22 .
84 ˚A per atom as a function of temperature for B, C, N and O. Forboth basal and c -axis diffusion, the activation volume of O is isotropic and negative below 740Kwhile it remains positive for B, C and N. The activation volume for all the solutes increases withincreasing temperatures, in part, as the elastic constants soften as temperature increases . Thisincrease is ∼
14% for basal activation volume and ∼
15% for c -axis activation volume for all thesolutes at 923K. where D ij (0) is the diffusivity tensor components at zero pressure. The activation volumeis calculated using V ij = − ( D ij (0)) − k B T ∂D ij ∂p (cid:12)(cid:12)(cid:12)(cid:12) p =0 = ( D ij (0)) − k B T (cid:88) kl d ijkk S kkll (17)where d is the elastodiffusivity tensor and S is the elastic compliance tensor. In the case ofinterstitial diffusion, the activation volume is equal to the migration volume of a jump: thevolume change between the transition state and initial state .Figure 6 shows that the activation volume for O diffusion is isotropic and negative below740K, which leads to an increase in basal and c -axis diffusivities under hydrostatic pressure.The activation volumes for B, C and N diffusion remain positive throughout the tempera-ture range, with N having the largest activation volume. For O diffusion below 740K, thedominating t-o transition has negative migration volume, while the dominating transitionsfor the diffusion of other solutes have positive migration volumes. Negative activation vol-19me has also been observed experimentally for C diffusing in hcp-Co and in α -Fe , andtheir magnitudes are comparable to the activation volume computed for O diffusion in Mg.Due to the temperature-induced softening of the elastic constants , the activation volumeof basal and c -axis diffusion increases by ∼
14% and ∼
15% from 300K to 923K for all foursolutes. /T . − 5 H O D W L Y H F K D Q J H L Q D ∆ D Bb ∆ D Cb ∆ D Nb ∆ D Ob ∆ D B c ∆ D C c ∆ D N c ∆ D O c T . FIG. 7. (color online) Change in basal and c -axis diffusivity due to thermal strain, relative to thestrain free diffusivity for B, C, N and O. The thermal strain is nearly isotropic and linear over theentire temperature range, to a maximum value of 2% at the melting temperature of 923K. Theeffect of thermal expansion is largest for N, for which the diffusivity doubles approaching melting,and smallest for O. Below 740K, O diffusivity decreases relative to its strain free diffusivity—dueto the negative activation volume—unlike the other three solutes. Thermal expansion effect on diffusion
Figure 7 shows that thermal expansion increases the diffusivity of B, C and N, butdecreases the diffusivity of O up to 740K. The fit of experimental thermal expansion datato temperature is used to estimate thermal strain. Thermal expansion is nearly isotropicin the temperature range 300K to 923K, reaching a maximum value of 2%. For B, C and N,both basal and c -axis diffusivities increase upon thermal expansion, with N experiencing thelargest effect of more than 100% increase in diffusivity at T > c -axis. Above 740K the O diffusivity is greatercompared to its strain free diffusivity as expected due to thermal expansion. However, below740K the O diffusivity is lower compared to its strain free diffusivity. This non-montonicbehavior of O diffusivity with thermal expansion is due to the sign inversion of five of theelastodiffusion tensor components. VI. CONCLUSION
We determine the stable interstitial sites, migration barriers, diffusivities, and elasto-diffusion tensors for B, C, N and O in Mg. We find a new stable distorted hexahedral sitethat B and C can occupy in Mg. Analytical expressions for interstitial diffusion in bulkhcp crystals are derived for the networks of interstitial sites. Diffusion of O is isotropic dueto dominating isotropic t-o transitions while B and C have faster basal diffusion comparedto c -axis diffusion and N have slower basal diffusion compared to c -axis diffusion. Thisshows that diffusion depends on the diffusion network formed by sites and their energetics,which varies from solute to solute. The elastodiffusion tensor captures the effect of strainon diffusivity by summing the contributions from changes in jump vectors and changes inmigration barriers. For B, C, N and O in Mg, the contribution to elastodiffusion compo-nents due to changes in migration barriers dominates over the contribution from changesin jump vectors with a few exceptions. There are a few elastodiffusion components whichchange sign at crossover temperature due to competing mechanisms. In the case of O, fiveof the elastodiffusion components change sign, which leads to negative activation volumebelow 740K and decreased diffusivity upon thermal expansion. This behavior of O as aninterstitial defect is counterintuitive because interstitial diffusivity is expected to decreaseunder compression as transition states are usually “smaller.” We see that N in its groundstate (octahedral) contracts the crystal upon relaxation while it has the positive activationvolume; O in its ground state (tetrahedral) expands the crystal on relaxation while hav-ing a negative activation volume. This shows that elastic dipole tensor of transition statesplays a vital role along with the energetics of sites. Our study of interstitial solute diffusionunder strain can be extended for other crystal structures and interstitial defects. Finally,understanding interstitial solute kinetics under strain can be helpful in studying the solutediffusivity in the heterogeneous strain fields due to dislocations or other defects.21 CKNOWLEDGMENTS
Figures 1 and 2 are generated using the Jmol package . The authors thank GraemeHenkelman for helpful conversations. This research was supported by the U.S. Office ofNaval Research under the grant N000141210752 and the National Science Foundation Award1411106, with computing resources provided by the University of Illinois campus clusterprogram. The full tabular data is archived by NIST at materialsdata.nist.gov ; seeRef. 26. The authors also thank the Library Service at Los Alamos National Laboratoryfor locating a copy of Ref. 10, and Yulia Maximenko at Univ. Illinois, Urbana-Champaign,Dept. of Physics for her translation help. ∗ [email protected] T. M. Pollock, Science , 986 (2010). H. E. Friedrich and B. L. Mordike,
Magnesium Technology , 1st ed. (Springer-Verlag BerlinHeidelberg, 2006). W. Joost, JOM , 1032 (2012). E. Fromm and G. H¨orz, International Metals Reviews , 269 (1980). X.-Z. Wu, L.-L. Liu, R. Wang, L.-Y. Gan, and Q. Liu, Frontiers of Materials Science , 405(2013). Y. Bertin, J. Parisot, and J. Gacougnolle, Journal of the Less Common Metals , 121 (1980). H. H. Wu and D. R. Trinkle, Phys. Rev. Lett. , 045504 (2011). A. O’Hara and A. A. Demkov, Applied Physics Letters , 211909 (2014). H. H. Wu, P. Wisesa, and D. R. Trinkle, Phys. Rev. B , 014307 (2016). V. Zotov and A. Tseldkin, Soviet Physics Journal , 1652 (1976). D. R. Trinkle, Philos. Mag. (2016), 10.1080/14786435.2016.1212175. D. R. Trinkle, “
Onsager ,” http://dallastrinkle.github.io/Onsager (2016). P. H. Dederichs and K. Schroeder, Phys. Rev. B , 2524 (1978). E. J. Savino and N. Smetniansky-De Grande, Phys. Rev. B , 6064 (1987). C. H. Woo and C. B. So, Philosophical Magazine A , 1299 (2000). G. Kresse and J. Furthm¨uller, Phys. Rev. B , 11169 (1996). P. E. Bl¨ochl, Phys. Rev. B , 17953 (1994). G. Kresse and D. Joubert, Phys. Rev. B , 1758 (1999). J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. , 3865 (1996). M. Methfessel and A. T. Paxton, Phys. Rev. B , 3616 (1989). J. Friis, G. K. H. Madsen, F. K. Larsen, B. Jiang, K. Marthinsen, and R. Holmestad, TheJournal of Chemical Physics , 11359 (2003). G. Henkelman, B. P. Uberuaga, and H. J´onsson, The Journal of Chemical Physics , 9901(2000). This approximation introduces at most a 40% error in the attempt frequencies; the error isestimated by comparing with a large Mg supercell using bulk force constants, and introducingthe interstitial-Mg force constants from the finite displacement calculations. G. H. Vineyard, Journal of Physics and Chemistry of Solids , 121 (1957). Following Varvenne et al. , we can estimate the finite-size error from using a 4 × × R. Agarwal and D. R. Trinkle, “Data citation: Light element diffusion in mg using first principlescalculations: Anisotropy and elastodiffusion,” (2016), http://hdl.handle.net/11256/694. E. Clouet, S. Garruchet, H. Nguyen, M. Perez, and C. S. Becquart, Acta Materialia , 3450(2008). H. Mehrer, Defect and Diffusion Forum , 91 (2011). C. W. Greeff and J. A. Moriarty, Phys. Rev. B , 3427 (1999). M. Wuttig and J. Keiser, Phys. Rev. B , 815 (1971). A. Bosman, P. Brommer, L. Eijkelenboom, C. Schinkel, and G. Rathenau, Physica , 533(1960). Y. Touloukian, R. Kirby, R. Taylor, and T. Lee,
Thermophysical Properties of Matter: ThermalExpansion-Metallic Elements and Alloys , 1st ed., Vol. 12 (Plenum Press, New York, 1975). Jmol: an open-source Java viewer for chemical structures in 3D. C. Varvenne, F. Bruneval, M.-C. Marinica, and E. Clouet, Phys. Rev. B , 134102 (2013)., 134102 (2013).