Diffusivity and derivatives for interstitial solutes: Activation energy, volume, and elastodiffusion tensors
JJuly 13, 2016 Philosophical Magazine paper
To appear in
Philosophical Magazine
Vol. 00, No. 00, Month 20XX, 1–22 Di ff usivity and derivatives for interstitial solutes: Activation energy,volume, and elastodi ff usion tensors Dallas R. Trinkle ∗ Department of Materials Science and Engineering, University of Illinois, Urbana-Champaign ( July 13, 2016 ) Computational atomic-scale methods continue to provide new information about geometry, ener-getics, and transition states for interstitial elements in crystalline lattices. This data can be used todetermine the di ff usivity of interstitials by finding steady-state solutions to the master equation.In addition, atomic-scale computations can provide not just the site energy, but also the stressin the cell due to the introduction of the defect to compute the elastic dipole. We derive a gen-eral expression for the fully anistropic di ff usivity tensor from site and transition state energies,and three derivatives of the di ff usivity: the elastodi ff usion tensor (derivative of di ff usivity withrespect to strain), the activation barrier tensor (logarithmic derivative of di ff usivity with respectto inverse temperature) and activation volume tensor (logarithmic derivative of di ff usivity withrespect to pressure). Computation of these quantities takes advantage of crystalline symmetry,and we provide an open-source implementation of the algorithm. We provide analytic resultsfor octahedral-tetrahedral networks in face-centered cubic, body-centered cubic, and hexagonalclosed-packed lattices, and conclude with numerical results for C in Fe. Keywords:
Interstitial di ff usion; activation barrier; elastodi ff usion tensor; automatedcomputation
1. Introduction
Mass transport in solids is an integral part of material processing for materials from metalsto ceramics to semiconductors. At the atomic scale, atoms move via defects: either va-cancies for atoms on the crystalline lattice, or via interstitial sites o ff of the lattice[1, 2].For the case of interstitial defects, what had been often considered the “simplest” trans-port to model can often conceal surprising complexity[3] when the interstitial networkinvolves more than a single symmetry unique site. Furthermore, real materials also havenon-homogeneous strain fields, and strain provides both a driving force[4] for di ff usionand modifies transport coe ffi cients leading to a variety of complex behavior including:anisotropic dopant di ff usion in semiconductor thin films[5–7], internal stress fields near adislocation a ff ecting transport and segregation of C solute atoms in Fe[8, 9] or vacanciesin HCP metals[10], and anisotropy of dumbbell interstitial di ff usion under biaxial stress inCu and Pt[11]. The influence of strain on di ff usivity has been previously investigated onlyfor cases where the interstitial random walk is purely uncorrelated[12, 13]; this excludeseven the simple case of octahedral-tetrahedral networks in hexagonal-closed packed mate-rials. While advances in density functional theory using modern supercomputers combinedwith transition-state finding methods[14, 15] can compute the energies for di ff erent con- ∗ Email: [email protected] a r X i v : . [ c ond - m a t . m t r l - s c i ] J u l uly 13, 2016 Philosophical Magazine paper figurations and transition rates at the atomic scale with near chemical accuracy[16, 17],general approaches to (a) automate the generation of a transition network, and (b) computetransport coe ffi cients and related quantities using non-stochastic techniques are lacking.The simplest transport coe ffi cient to consider is the (generally anisotropic) di ff usivitytensor D , but we can also conceive of derivatives with respect to temperature and strain (orstress). An interstitial solute executes a random walk over a connected network of sites,where the transition rate between sites is given by a rate following an Arrhenius relationfrom transition-state theory[18]. From the long-time behavior, the di ff usivity tensor D canbe computed, where D relates a steady-state flux j to a concentration gradient −∇ c . Theexistence of di ff erent symmetry unrelated sites (and transitions) means that the total di ff u-sivity may not have a simple Arrhenius temperature dependence given by a single activa-tion barrier. Instead, over a su ffi ciently small temperature range, if we were to approximatethe di ff usivity tensor with an Arrhenius relation, D ≈ D exp( − β E act ), then the negativederivative of di ff usivity with respect to inverse temperature β allows us to define E act : = − D − / dDd β D − / , (1)corresponding to an anisotropic (and possibly temperature dependent) “barrier” tensor.The barrier should correspond to the activation energy for the rate-limiting transition indi ff usion. The fourth-rank elastodi ff usion tensor[1, 13], d : = dDd ε , (2)is the derivative of di ff usivity with respect to strain. As both di ff usivity and strain aresymmetric tensors, d abcd = d bacd = d abdc ; furthermore, d obeys crystalline symmetry.Thus, it has a similar structure to elastic sti ff nesses and compliances, with one excep-tion: as it is not the second derivative of an energy, there is no general relationship be-tween d abcd = dD ab / d ε cd and d cdab = dD cd / d ε ab , unless imposed by crystalline symmetry.Finally, activation volume can be defined in a similar fashion to activation energy, as alogarithmic derivative with respect to stress, V act abcd : = k B T D − / dDd σ D − / = k B T (cid:88) i jkl = (cid:16) D − / (cid:17) ai dD i j d ε kl (cid:16) D − / (cid:17) jb S klcd , (3)and related back to the elastodi ff usion tensor using the elastic compliances S i jkl . Thisfourth-rank tensor contains the full aniostropy of di ff usion response to stress. It can bereduced to a scalar activation volume by a double contraction, V act : = (cid:88) i j = V act ii j j (4)where the second contraction considers hydrostatic pressure, and the first contraction aver-ages in all spatial directions.In what follows, we develop the theory for the di ff usivity, elastodi ff usivity, activationbarrier and volume tensors of an interstitial in terms of the site energies, transition stateenergies, and elastic dipoles (derivative of energy with respect to strain) of those same2 uly 13, 2016 Philosophical Magazine paper quantities. This includes the cases where the interstitial random walk includes correla-tion. We apply the theory to three simple examples: octahedral-tetrahedral networks inface-centered cubic, body-centered cubic, and hexagonal closed-packed lattices. We alsoprovide numerical values for carbon di ff usion in iron, using data from[9]. Finally, the workoutlined here for the general case is implemented in open-source, publicly-available code(c.f. Appendix C and [19]) that can be applied to other systems[20].
2. Methodology2.1. Di ff usivity In the case of interstitial di ff usivity in the dilute limit, an interstitial atom moves through anetwork of sites in the crystal. Such a model allows simplification of the master equation:the system state is full described by the location of the interstitial, x = R + u for R a latticevector and u a vector in the unit cell. Due to translational invariance, the equilibrium siteprobability and set of transitions from a site are determined by the vector in the unit cell.For a site i in the unit cell, it has an equilibrium site probability ρ i that follows an Arrheniusrelationship, ρ i : = Z − ρ i exp ( − β E i ) (5)for site energy E i , entropic prefactor ρ i = exp( S i / k B ) , and partition function Z = (cid:80) i ρ i exp ( − β E i ). The transition from site i to site j has a rate λ i → j , λ i → j : = λ i j ρ i exp (cid:16) − β (cid:104) E ts i j − E i (cid:105)(cid:17) (6)for transition state energy E ts i j and entropic prefactor λ i j = exp( S ts i j / k B ), following [18]. Inthis formulation, the transition state energy and entropic prefactors are equal for i → j andfor j → i , while it is not necessary that λ i → j and λ j → i are equal. Finally, the probabili-ties obey detailed balance, where ρ i λ i → j = ρ j λ j → i for all i , j . Each transition results in adisplacement of the di ff using atom by δ x i → j : = x j − x i . (7)With these definitions of site probabilities and rates, we can write down an expression forthe di ff usivity[2]. There are two contributions to the di ff usivity: an uncorrelated di ff usivityand a correlation correction due to the unbalanced hops from individual sites. We definetwo rate matrices, Λ i j , Λ i j : = λ i → j : i (cid:44) j − (cid:80) j λ i → j : i = j (8)which is the transition matrix for the master equation, and its symmetric counterpart ω i j : = ρ / i Λ i j ρ − / j . (9)These matrices are negative-definite, and—due to detailed balance—both have a null-spacerelated to the site probabilities ρ i . We define the scaled velocity vector (corresponding to3 uly 13, 2016 Philosophical Magazine paper the bias of jumps at a site as in [8, 9]), b i : = ρ / i (cid:88) j λ i → j δ x i → j (10)which is non-zero when there are unbalanced hops from a site i ; together with the bias-correction vector γ i which solves (cid:88) j ω i j γ j = b i . (11)This is most easily done with the pseudoinverse of ω , g , so that γ i = (cid:80) j g i j b j . Finally, thedi ff usivity is D = (cid:88) i j δ x i → j ⊗ δ x i → j λ i → j ρ i + (cid:88) i b i ⊗ γ i (12)where ⊗ is the outer (or dyad) product of two vectors. Because ω is a symmetric matrix,so is g , and hence D is a symmetric second-rank tensor.While the expressions leading up to Eqn. 12 are general for non-interacting atoms hop-ping between sites, crystal symmetry can reduce the complexity of the expressions to beevaluated. We use the Seitz notation[21] for a symmetry operation { R , t } , where for a point x , { R , t } x : = R x + t . Then, the inverse { R , t } − = { R − , − R − t } . The full set of operationsmake up the space group , and at any site i there will be the subgroup of operations thatleave that site fixed, its point group . If every site i has a position x i , then each group oper-ation can be expressed as a matrix { R , t } i j : = δ ( x i − { R , t } x j ) (13)where δ is the Kronecker delta. This matrix then applies the group operation { R , t } to a sitescalar f i as (cid:80) j { R , t } i j f j . A related matrix can be defined corresponding to vector at eachsite (such as b i and γ i ), where with the orthonormal basis vectors { e a } , then the matrix { R , t } ia , jb : = { R , t } i j ( e a · R · e b ) (14)transforms any site vector f i as (cid:80) a e a (cid:80) jb { R , t } ia , jb ( f i · e b ). All of our site scalars ( ρ i ), sitevectors ( b i and γ i ), and site-to-site matrices ( Λ i j , ω i j , and g i j ) share the crystal symmetry, The construction a ⊗ b is a second rank tensor such that ( a ⊗ b ) · v = ( b · v ) a for any vector v . uly 13, 2016 Philosophical Magazine paper and so the following symmetry relations are true for any space group operation { R , t } , ρ i = (cid:88) j { R , t } i j ρ j , b i = (cid:88) a e a (cid:88) jb { R , t } ia , jb ( b j · e b ) , γ i = (cid:88) a e a (cid:88) jb { R , t } ia , jb ( γ j · e b ) , Λ i j = (cid:88) kl { R , t } ik Λ kl { R , t } l j ,ω i j = (cid:88) kl { R , t } ik ω kl { R , t } l j , g i j = (cid:88) kl { R , t } ik g kl { R , t } l j . (15)The implication of these symmetries is that ρ i , b i , and γ i can be written as linear combina-tions of vectors left unchanged by all { R , t } matrices, and that our matrices Λ , ω , and g canbe entirely expanded in that same basis.The basis functions for site scalars and site vectors come from the Wycko ff positions cor-responding to sites in the network. Each Wycko ff position represents a full set of symmetry-related sites; any site scalar that obeys crystal symmetry will have the same value for eachsite corresponding to the same Wycko ff position. Hence, if we have a Wycko ff positionthat has N W sites in the unit cell, its basis vector components will be N − / for each sitecorresponding to that Wycko ff position, and 0 otherwise; there will be one basis func-tion for each unique Wycko ff site regardless of how many sites that may represent in theunit cell . The Wycko ff positions also serve to construct a basis for site vectors combinedwith the point group operations for that site. Considering the point group operations thatare available in a crystal, there are ten types: identity, a 2-, 3-, 4-, or 6-fold axis, mirrorthrough a plane, or a 2-, 3-, 4-, or 6-fold axis combined with a mirror through that sameplane. Note that a 2-fold axis combined with a mirror operation is inversion. As we areinterested in vectors that remain unchanged by all of the point group operations for a site,we can easily generate the basis corresponding to each operation. For identity, we have any3-dimensional basis; for a mirror operation, a 2-dimensional basis that spans the mirrorplane; for a 2-, 3-, 4-, or 6-fold axis, a 1-dimensional basis corresponding to the axis; andfor a 2-, 3-, 4-, or 6-fold axis combined with a mirror, there are no vectors left unchanged.We construct the vector space left unchanged by each operation, and intersect to generatethe final vector space for a single site. The spanning vectors for that site can then be rotatedto the other sites corresponding to the same Wycko ff position, and normalized in a similarway to the site scalar basis.The crystal symmetry divides di ff usion networks into three types based on the crystalsymmetry: networks without correlation, networks with correlation and inversion, and net-works with correlation but without inversion. As it is possible for the site vector basis tobe an empty set (e.g., if each Wycko ff position point group possesses a 2-, 3-, 4-, or 6-foldaxis combined with a mirror operation), this excludes correlation explicitly by symmetry .Next, there may be a non-empty site vector basis while the crystal has inversion. In thiscase, the transition matrices Λ or ω expressed in the site vector basis are not singular. Thisis because the null-space vector of Λ —which corresponds to ρ i —has zero projection intoany basis vector. Finally, in the most general case, the transition matrices remain singular,5 uly 13, 2016 Philosophical Magazine paper and so the construction of g requires the pseudoinverse of ω . Activation barrier
We define the activation energy tensor from the first derivative of the di ff usivity with re-spect to inverse temperature, and evaluate with perturbation theory. In this case, we intro-duce a small change in inverse temperature, − δβ , to ρ i and λ i → j , then evaluate the first orderchange in D . Then, the first order change in ρ i , δρ i , is given by ρ i + δρ i = ρ i e δβ E i (cid:80) j ρ j e δβ E j = ρ i + ρ i ( E i − (cid:104) E (cid:105) ) δβ + O ( δβ ) (16)where the average energy (cid:104) E (cid:105) = (cid:80) i ρ i E i . Next, the first order change in λ i → j , δλ i → j , isgiven by λ i → j + δλ i → j = λ i → j e δβ (cid:16) E ts ij − E i (cid:17) = λ i → j + λ i → j (cid:16) E ts i j − E i (cid:17) δβ + O ( δβ ) . (17)From these, we can compute the related terms: δ b i and δω . First, b i + δ b i = ( ρ i + δρ i ) / (cid:88) j ( λ i → j + δλ i → j ) δ x i → j = b i + ρ / i (cid:88) j λ i → j (cid:34) E ts i j −
12 ( E i + (cid:104) E (cid:105) ) (cid:35) δ x i → j δβ + O ( δβ ) . (18)Next for δω ii we have δω ii = − (cid:88) j δλ i → j = − (cid:88) j λ i → j (cid:16) E ts i j − E i (cid:17) δβ + O ( δβ ) (19)and for δω i j ( i (cid:44) j ) we have δω i j = ( ρ i + δρ i ) / ( λ i → j + δλ i → j )( ρ j + δρ j ) − / = ω i j (cid:34) E ts i j − (cid:16) E i + E j (cid:17)(cid:35) δβ + O ( δβ ) . (20)Finally, the first-order change in the pseudoinverse g , δ g is given by (c.f. Appendix A) g + δ g = ( ω + δω ) − = g − g δω g + O ( δβ ) . (21)6 uly 13, 2016 Philosophical Magazine paper With all of these terms, we can compute the first-order correction to the di ff usivity, whichgives the derivative of D with respect to − β by substituting into Eqn. 12, − dDd β = (cid:88) i j δ x i → j ⊗ δ x i → j (cid:104) E ts i j − (cid:104) E (cid:105) (cid:105) λ i → j ρ i + (cid:88) i (cid:32) δ b i δβ ⊗ γ i + γ i ⊗ δ b i δβ (cid:33) − (cid:88) i j γ i ⊗ δω i j δβ γ j (22)where the first term is the uncorrelated contribution and the remaining terms come fromcorrelation. It should be noted that all of the perturbed site scalar, site vector, and matriceshave the same symmetry as their unperturbed versions. A similarly structured solution willfollow for the evaluation of the elastodi ff usion tensor. Elastodi ff usion tensor To evaluate the elastodi ff usion tensor, we introduce a perturbation through strain whichchanges site and transition state energies, as defined by the elastic dipole tensor. The elasticdipole tensor P i for a site i is P i : = − dE i d ε . (23)The elastic dipole can be conveniently evaluated in a supercell calculation from the stressin the cell: an interstitial is added to an initially undefected, unstressed supercell containing N atoms (with equilibrium volume V per atom), resulting in a stress σ , then to first orderin N − , P ≈ NV σ, (24)which is straightforward to evaluate with density-functional theory methods; e.g., see [22–25]. Similarly, the energy of a transition state can also change with strain, as dictated bythe elastic dipole tensor for the transition state P ts i j for the transition state between i and j , P ts i j : = − dE ts i j d ε . (25)This, too, can be approximated by the stress at the transition state in a supercell calculationas in Eqn. 24; e.g., see [9, 24]. The definitions of elastic dipoles allow the introduction ofa small strain perturbation δε to produce site energies changes δ E i and transition energies δ E ts i j as δ E i = − P i : δεδ E ts i j = − P ts i j : δε (26)which is correct to first order in strain. The double contraction sums over both indices of the two second-rank tensors: A : B = (cid:80) ab A ab B ba . uly 13, 2016 Philosophical Magazine paper We define the elastodi ff usion tensor from the first derivative of the di ff usivity with re-spect to strain, and evaluate with perturbation theory. Strain produces changes in site en-ergies and transition state energies as laid out above, but also modifies all vectors in thenetwork. The purely geometric contribution is evaluated by expanding the di ff usivity tensorout in components as D = (cid:80) ab D ab e a ⊗ e b for an orthonormal basis { e a } ; then a (symmetric)strain perturbs the basis vectors by δ e a = (cid:80) c δε ac e c , so that δ D (geom) : = (cid:88) ab e a ⊗ e b (cid:88) c D ac δε cb + D bc δε ca (27)which produces a (symmetrized) contribution to the elastodi ff usion tensor d (geom) abcd : =
12 ( δ ad D bc + δ ac D bd + δ bc D ad + δ bd D ac ) . (28)To get the uncorrelated and correlated changes that we add to this, we get the first ordercorrection to ρ i , δρ i as ρ i + δρ i = ρ i e − βδ E i (cid:80) j ρ j e − βδ E j = ρ i + ρ i β (cid:16) P i − (cid:104) P (cid:105) (cid:17) : δε + O ( δε ) (29)where the average elastic dipole (cid:104) P (cid:105) = (cid:80) i ρ i P i . Next, the first order change in λ i → j , δλ i → j ,is given by λ i → j + δλ i → j = λ i → j e − β (cid:16) δ E ts ij − δ E i (cid:17) = λ i → j + λ i → j β (cid:16) P ts i j − P i (cid:17) : δε + O ( δε ) . (30)From these, we can compute the related terms: δ b i —ignoring contributions of the strain tothe hop vectors—and δω . First, b i + δ b i = ( ρ i + δρ i ) / (cid:88) j ( λ i → j + δλ i → j ) δ x i → j = b i + ρ / i (cid:88) j δ x i → j λ i → j β (cid:34) P ts i j − (cid:16) P i + (cid:104) P (cid:105) (cid:17)(cid:35) : δε + O ( δε ) . (31)Next for δω ii we have δω ii = − (cid:88) j δλ i → j = − (cid:88) j λ i → j β (cid:16) P ts i j − P i (cid:17) : δε + O ( δε ) (32)and for δω i j ( i (cid:44) j ) we have δω i j = ( ρ i + δρ i ) / ( λ i → j + δλ i → j )( ρ j + δρ j ) − / = ω i j β (cid:34) P ts i j − (cid:16) P i + P j (cid:17)(cid:35) : δε + O ( δε ) . (33)8 uly 13, 2016 Philosophical Magazine paper The expression for the first-order change in the pseudoinverse is the same as Eqn. 21;hence, we can combine all of the contributions into Eqn. 12 (similar to Eqn. 22) plusEqn. 28 to get d abcd =
12 ( δ ad D bc + δ ac D bd + δ bc D ad + δ bd D ac ) + (cid:88) i j δ x i → j , a δ x i → j , b β (cid:16) P ts i j , cd − (cid:104) P (cid:105) cd (cid:17) λ i → j ρ i + (cid:88) i (cid:32) δ b i , a δε cd γ i , b + γ i , a δ b i , b δε cd (cid:33) − (cid:88) i j γ i , a δω i j δε cd γ j , b (34)where the first term is the geometric contribution, the second term is the uncorrelated con-tribution, and the remaining terms come from correlation.Introducing a finite strain into the lattice can reduce the symmetry of the lattice, and cou-ples with the anisotropy of the elastic dipole tensors. The dipole tensors can be expandedin a site tensor basis, similar to the site vector basis, using the point group for a singlerepresentative site of each Wycko ff position. We can also make a similar expansion for thetransition states, by working with the “double point group”: the set of symmetry operationsthat leave both the initial and final positions unchanged. Despite the breaking of symmetry,we can see from Eqn. 34, that if a site i has an empty vector basis (so that both b i and γ i are zero by symmetry), the correlation contributions to the elastodi ff usion tensor are also zero. This means that the introduction of strain does not require the computation of δ b i or δω i j for any sites where b i and ω i j was not already computed. However, it is possible thatparticular values of δε can produce δ b i with components that would be zero by symmetryfor b i ; the symmetry of δ b i /δε jk is that of a third-rank tensor at site i .
3. Results
We can apply these results to three classic systems: octahedral-tetrahedral networks in face-centered cubic (FCC), body-centered cubic (BCC), and hexagonal closed-packed (HCP)structures. In the case of FCC and BCC, the network has zero site vectors basis (no cor-relation contribution), while the HCP network has a site vector basis oriented along the c -axis. In the case of BCC, tetrahedrals can be transition states, or states on their own.These three cases move from one, to two, to three di ff erent transition state energies. Fi-nally, we conclude with a short set of numerical results for carbon in iron, using the EAMdata generated by [9]. Face-centered cubic
In a face-centered cubic lattice, the common interstitial sites are octahedral sites that candi ff use to tetrahedral sites, and vice versa. The face-centered cubic space group is Fm ¯3 m ,where the solvent atoms occupy 4 a Wycko ff positions, with octahedrals at 4 b and tetra-hedrals at 8 c .[26] There is one octahedral site for every “solvent” atom, forming an inter-penetrating FCC lattice. Each octahedral is connected to eight tetrahedral sites via a (cid:104) (cid:105) jumps. There are two tetrahedral sites for every lattice atom, forming an interpenetratingsimple cubic lattice. Each tetrahedral is connected to four octahedral sites; the jump vectorsare a subset of a (cid:104) (cid:105) vectors with tetrahedral symmetry (either { [111] , [1¯1¯1] , [¯11¯1] , [¯1¯11] } or { [¯1¯1¯1] , [¯111] , [1¯11] , [11¯1] } ). By symmetry, there is a single transition state energy, E tsot ,9 uly 13, 2016 Philosophical Magazine paper so that we have two rates λ o → t ∝ exp (cid:104) − β (cid:16) E tsot − E o (cid:17)(cid:105) λ t → o ∝ exp (cid:104) − β (cid:16) E tsot − E t (cid:17)(cid:105) (35)and we have the probabilities ρ o = λ t → o λ o → t + λ t → o , ρ t = λ o → t λ o → t + λ t → o , (36)such that ρ o + ρ t =
1. The site vector basis for all of our sites are empty, so that we haveno correlation terms to consider. From this, we have the isotropic di ff usivity from Eqn. 12, D = (cid:40) a λ o → t ρ o + · a λ t → o ρ t (cid:41) = a λ o → t ρ o = a λ t → o ρ t , (37)and the isotropic activation energy from Eqn. 1 and Eqn. 22, E act = D − (cid:40) a (cid:104) E tsot − (cid:104) E (cid:105) (cid:105) λ o → t ρ o + · a (cid:104) E tsot − (cid:104) E (cid:105) (cid:105) λ t → o ρ t (cid:41) = (cid:104) E tsot − (cid:104) E (cid:105) (cid:105) . (38)We can also determine the elastodi ff usion tensor and (isotropic) activation volume foran interstitial in FCC. To determine the elastodi ff usion tensor, we note that the octahedraland tetrahedral elastic dipoles are fully isotropic by cubic and tetrahedral symmetry, re-spectively. The elastic dipole corresponding to the transition state has a three-fold rotationaxis given by the (cid:104) (cid:105) vector connecting the two sites, so that the dipole has a componentparallel P tsot , (cid:107) to the rotation axis v = √ (cid:104) (cid:105) and a component perpendicular P tsot , ⊥ , P tsot ( v ) = P tsot , (cid:107) v ⊗ v + P tsot , ⊥ ( − v ⊗ v ) = P tsot , ⊥ + ( P tsot , (cid:107) − P tsot , ⊥ ) v ⊗ v (39)which, in component form, is P tsot ([ s , s , s ]) ab = P tsot , ⊥ δ ab + (cid:16) P tsot , (cid:107) − P tsot , ⊥ (cid:17) s a s b s a = ± v . Similarly, ( δ x ⊗ δ x ) ab = s a s b a /
16. We can then evaluate the elastodi ff usion tensor, using the Voigt notation forfourth-rank tensors, and Eqn. 34, d = d = d = a (cid:40) + β (cid:32) P tsot , ⊥ + P tsot , (cid:107) − (cid:104) P (cid:105) (cid:33)(cid:41) λ o → t ρ o d = d = d = d = d = d = a (cid:40) β (cid:32) P tsot , ⊥ + P tsot , (cid:107) − (cid:104) P (cid:105) (cid:33)(cid:41) λ o → t ρ o d = d = d = a (cid:40) + β (cid:32) P tsot , (cid:107) − P tsot , ⊥ (cid:33)(cid:41) λ o → t ρ o (41)10 uly 13, 2016 Philosophical Magazine paper and all other terms are zero, and where (cid:104) P (cid:105) = P o ρ o + P t ρ t . Due to symmetry, both P o and P t are isotropic (scalars). Then, the activation volume tensor from Eqn. 3 is V act11 = k B T S + (cid:32)
13 Tr P tsot − (cid:104) P (cid:105) (cid:33) B − V act12 = k B T S + (cid:32)
13 Tr P tsot − (cid:104) P (cid:105) (cid:33) B − V act44 = k B T S +
16 ( P tsot , (cid:107) − P tsot , ⊥ ) S (42)for bulk modulus B − = S + S ), and where Tr P tsot = P tsot , ⊥ + P tsot , (cid:107) . The terms propor-tional to k B T all represent geometric contributions, while the remaining terms correspondto the sensitivity of the transition states and interstitial sites to stress. Body-centered cubic
In a body-centered cubic crystal, the common interstitial site is an octahedral site, withthe possibility of tetrahedral sites either as metastable or transition sites. Here, to com-pare with the face-centered cubic and hexagonal closed-packed case, we will take bothoctahedral and tetrahedral sites to be stable, and so there will be di ff usion from octahedralto tetrahedral sites and between tetrahedral sites. The body-centered cubic space group is Im ¯3 m , where the solvent atoms occupy 2 a Wycko ff positions, with octahedrals at 6 b andtetrahedrals at 12 d .[26] Unlike our other cases, there are three octahedral sites and sixtetrahedral sites for every “solvent” atom; moreover, the octahedral sites do not have cu-bic point group symmetry, but rather tetragonal symmetry. If we consider a representativeoctahedral at [00 ], it connects to four tetrahedrals in the (001) plane, with jump vectors ± a [100] and ± a [010] and a transition state energy E tsot . The elastic dipole associated withthis octahedral has two components, P o = P o , ⊥ ( x ⊗ x + y ⊗ y ) + P o , p z ⊗ z , corresponding to the in-plane P o , p and perpendicular P o , ⊥ components. The transition statealong [100] has an elastic dipole with three components, P tsot = P tsot , (cid:107) x ⊗ x + P tsot , ⊥ y ⊗ y + P tsot , p z ⊗ z , corresponding to the along-hop P tsot , (cid:107) as well as perpendicular and in-plane components.Considering next a representative tetrahedral at [0
14 12 ], it connects to two octahedrals andfour tetrahedrals, with jump vectors (cid:110) a [01¯1] , a [011] , a [¯110] , a [1¯10] (cid:111) and a transition stateenergy E tstt . The elastic dipole associated with this tetrahedral has two components, P t = P t , ⊥ ( x ⊗ x + z ⊗ z ) + P t , a y ⊗ y , corresponding to the four-fold rotation / mirror axis P t , a and perpendicular P t , ⊥ components.The transition state along [01¯1] has an elastic dipole with three components, P tstt = P tstt , (cid:107)
12 ( y − z ) ⊗ ( y − z ) + P tstt , ⊥
12 ( y + z ) ⊗ ( y + z ) + P tstt , p x ⊗ x , uly 13, 2016 Philosophical Magazine paper corresponding to the along-hop as well as perpendicular and in-plane components. In orderfor the site probabilities to match what was found in FCC, Eqn. 36, we need to account forthe increased multiplicity of octahedral and tetrahedral sites. We choose to keep the samenormalization ρ o + ρ t =
1, and include a factor of 1 / ff ects in the computation of the di ff usivity or its derivatives.From this, we have the isotropic di ff usivity of Eqn. 12, D = a { ρ o λ o → t + ρ t λ t → t } , (43)and the isotropic activation energy from Eqn. 1 and Eqn. 22, E act = D − (cid:40) a (cid:104) E tsot − (cid:104) E (cid:105) (cid:105) ρ o λ o → t + a (cid:104) E tstt − (cid:104) E (cid:105) (cid:105) ρ t λ t → t (cid:41) = (cid:2) E tsot − (cid:104) E (cid:105) (cid:3) ρ o λ o → t + (cid:2) E tstt − (cid:104) E (cid:105) (cid:3) ρ t λ t → t ρ o λ o → t + ρ t λ t → t . (44)The activation energy changes from E tsot −(cid:104) E (cid:105) when the octahedral to tetrahedral rate ρ o λ o → t is the fastest to E tstt −(cid:104) E (cid:105) when the tetrahedral to tetrahedral rate ρ t λ t → t is the fastest. Finally,the elastodi ff usion tensor has the cubic symmetry of the lattice, and using the Voigt notationfor fourth-rank tensors with Eqn. 34, d = d = d = a (cid:26) ρ o λ o → t + ρ t λ t → t ) + β (cid:16) P tsot , (cid:107) − (cid:104) P (cid:105) (cid:17) ρ o λ o → t + β (cid:16) P tstt , (cid:107) + P tstt , ⊥ − (cid:104) P (cid:105) (cid:17) ρ t λ t → t (cid:27) d = d = d = d = d = d = a (cid:40) β (cid:16) P tsot , ⊥ + P tsot , p − (cid:104) P (cid:105) (cid:17) ρ o λ o → t + β (cid:32) P tstt , (cid:107) + P tstt , ⊥ + P tstt , p − (cid:104) P (cid:105) (cid:33) ρ t λ t → t (cid:41) d = d = d = a (cid:40) ( ρ o λ o → t + ρ t λ t → t ) + β (cid:32) P tstt , (cid:107) − P tstt , ⊥ (cid:33) ρ t λ t → t (cid:41) (45)and all other terms are zero, and where (cid:104) P (cid:105) = (cid:110) Tr P o ρ o + Tr P t ρ t (cid:111) = (cid:40)(cid:32) P o , ⊥ + P o , p (cid:33) ρ o + (cid:32) P t , ⊥ + P t , a (cid:33) ρ t (cid:41) . (46)Due to symmetry, (cid:104) P (cid:105) is isotropic even though P o and P t are not. Then, the activation12 uly 13, 2016 Philosophical Magazine paper volume tensor from Eqn. 3 is V act11 = k B T S + (cid:110) P tsot , (cid:107) S + (cid:16) P tsot , ⊥ + P tsot , p (cid:17) S − (cid:104) P (cid:105) B − (cid:111) ρ o λ o → t ρ o λ o → t + ρ t λ t → t + (cid:40) (cid:16) P tstt , (cid:107) + P tstt , ⊥ (cid:17) ( S + S ) + P tstt , p S − (cid:104) P (cid:105) B − (cid:41) ρ t λ t → t ρ o λ o → t + ρ t λ t → t V act12 = k B T S + (cid:40)(cid:16) P tsot , (cid:107) + P tsot , ⊥ + P tsot , p (cid:17) S + (cid:32) P tsot , (cid:107) +
12 ( P tsot , ⊥ + P tsot , p ) (cid:33) S − (cid:104) P (cid:105) B − (cid:41) ρ o λ o → t ρ o λ o → t + ρ t λ t → t + (cid:40)(cid:16) P tstt , (cid:107) + P tstt , ⊥ + P tstt , p (cid:17) S + (cid:32) P tstt , (cid:107) + P tstt , ⊥ + P tstt , p (cid:33) S − (cid:104) P (cid:105) B − (cid:41) ρ t λ t → t ρ o λ o → t + ρ t λ t → t V act44 = k B T S + (cid:16) P tstt , (cid:107) − P tstt , ⊥ (cid:17) ρ t λ t → t ρ o λ o → t + ρ t λ t → t S (47)for bulk modulus B − = S + S ). The terms proportional to k B T all represent geomet-ric contributions, while the remaining terms correspond to the sensitivity of the transitionstates and interstitial sites to stress. There is more anisotropy in the activation volume forBCC compared with FCC, due to the lowered symmetry in the transition states. Hexagonal closed-packed
In a hexagonal closed-packed crystal, there are a variety of interstitial sites available[3, 20, 27–29], but to compare with the our cubic cases, we will consider an octahedral-tetrahedral network. The hexagonal closed-packed space group is P mmc , where the sol-vent atoms occupy 2 c Wycko ff positions, with octahedrals at 2 a and tetrahedrals at 4 f .[26]Similar to FCC, there is one octahedral and two tetrahedral sites for every “solvent” atom,but the connectivity is di ff erent. The octahedral sites form an interpenetrated simple hexag-onal lattice, where each octahedral connects to two octahedrals directly above and be-low ( ± c [0001] jumps) and six tetrahedrals (either {± a [1¯10 c a ] , ± a [01¯1 c a ] , ± a [¯101 c a ] } or {± a [¯110 c a ] , ± a [0¯11 c a ] , ± a [10¯1 c a ] } ). The tetrahedrals are each connected to one othertetrahedral (either c [0001] or c [000¯1]) and three tetrahedrals (opposing the octahedraljumps, and with the c -axis contribution negating the tetrahedral jump contribution). Thisintroduces three transition state energies, E tsot , E tstt , and E tsoo with corresponding rates amongoctahedrals and tetrahedrals. The site probabilities match what was found in FCC, Eqn. 36.From this, we can generate the anisotropic di ff usivity tensor (take x and y to be orthonor-mal basis vectors in the basal plane and z along the c -axis), D = c z ⊗ z ρ o λ o → o + (cid:32) a x ⊗ x + a y ⊗ y (cid:33) ρ o λ o → t + · (cid:32) c z ⊗ z (cid:33) ρ o λ o → t + c z ⊗ z ρ t λ t → t + b t ⊗ γ t (48)once we solve for the velocity and bias-correction vectors. The octahedral sites possessa threefold axis (pointing along z ) and a mirror plane perpendicular to that axis, so there The tetrahedral sites are not required to sit at a c / ff parameter z tocharacter this position in the 4 f site. However, in the case of di ff usion, only the long-range contribution matters, so thatchange in jump vectors due to z exactly cancels out, and we use z = / uly 13, 2016 Philosophical Magazine paper is no vector basis at those sites; the tetrahedral sites also possess the same threefold axis,but no mirror plane, and so a vector can point along z at the tetrahedral sites. A mirroroperation through the basal plane maps one tetrahedral sites to the other connected bya tetrahedral-to-tetrahedral jump; thus, the single site vector basis alternates ± z for eachsuccessive (basal) plane of tetrahedrals. In this basis, the symmetrized rate matrix only hasa single entry based on the escape rate from a tetrahedral and the tetrahedral-to-tetrahedraljump, ω = (cid:16) √ − √ (cid:17) · (cid:32) − λ t → t − λ t → o λ t → t λ t → t − λ t → t − λ t → o (cid:33) · √ − √ = − (2 λ t → t + λ t → o ) (49)and the velocity vector at a tetrahedral is ± (cid:16) λ t → o c − λ t → t c (cid:17) ρ / z . After projecting intothe site vector basis, and solving for the bias-correction vector, we have b t ⊗ γ t = − c
32 (3 λ t → o − λ t → t ) λ t → o + λ t → t ρ t z ⊗ z (50)and thus the di ff usivity from Eqn. 48 is D = a ( x ⊗ x + y ⊗ y ) ρ o λ o → t + c z ⊗ z ρ o λ o → o + c z ⊗ z ρ t λ t → o λ t → t λ t → o + λ t → t (51)Computing the activation energy tensor is complicated due to the correlation terms. Theuncorrelated contribution to inverse temperature derivative is − dDd β = c z ⊗ z ρ o λ o → o (cid:16) E tsoo − (cid:104) E (cid:105) (cid:17) + (cid:32) a x ⊗ x + a y ⊗ y (cid:33) ρ o λ o → t (cid:16) E tsot − (cid:104) E (cid:105) (cid:17) + · (cid:32) c z ⊗ z (cid:33) ρ o λ o → t (cid:16) E tsot − (cid:104) E (cid:105) (cid:17) + c z ⊗ z ρ t λ t → t (cid:16) E tstt − (cid:104) E (cid:105) (cid:17) (52)The change in the velocity vector for the tetrahedral site is δ b t = ± (cid:32) c λ t → o (cid:34) E tsot −
12 ( E t + (cid:104) E (cid:105) ) (cid:35) − c λ t → t (cid:34) E tstt −
12 ( E t + (cid:104) E (cid:105) ) (cid:35)(cid:33) ρ / z (53)and for the symmetric rate matrix δω = − (cid:16) λ t → t (cid:104) E tstt − E t (cid:105) + λ t → o (cid:104) E tsot − E t (cid:105)(cid:17) . (54)After we combine all of the terms into Eqn. 22, we have an activation energy tensor E act = (cid:16) E tsot − (cid:104) E (cid:105) (cid:17) ( x ⊗ x + y ⊗ y ) + (cid:32) ρ o λ o → o + ρ t λ t → o λ t → t λ t → o + λ t → t (cid:33) − (cid:34) (cid:110) E tsoo − (cid:104) E (cid:105) (cid:111) ρ o λ o → o + (cid:40) λ t → o λ t → o + λ t → t (cid:16) E tstt − (cid:104) E (cid:105) (cid:17) + λ t → t λ t → o + λ t → t (cid:16) E tsto − (cid:104) E (cid:105) (cid:17)(cid:41) ρ t λ t → o λ t → t λ t → o + λ t → t (cid:35) ( z ⊗ z ) . (55)14 uly 13, 2016 Philosophical Magazine paper For basal plane di ff usivity, the activation barrier is ( E tsot − (cid:104) E (cid:105) ), as expected. For di ff usionalong the c -axis, the expression in Eqn. 55 simplifies in several key limits. If ρ o λ o → o is thefastest, then E act ≈ ( E tsoo − (cid:104) E (cid:105) ); di ff usion is dominated by direct octahedral-to-octahedraljumps. If it is not, then the octahedral-tetrahedral network dominates; if λ t → t is rate-limiting(slower than λ t → o ), E act ≈ ( E tsot − (cid:104) E (cid:105) ), otherwise E act ≈ ( E tstt − (cid:104) E (cid:105) ). These important limitsare all captured in a single expression for the activation barrier, and can be computedautomatically.The elastic dipole terms break into basal and c -axis components, with a tetragonal com-ponent for the octahedral-tetrahedral jump. The elastic dipole associated with this octahe-dral has two components, due to hexagonal symmetry, P o = P o , b ( x ⊗ x + y ⊗ y ) + P o , c z ⊗ z , corresponding to the basal P o , b and c -axis P o , c components. The octahedral-to-octahedraltransition along the c -axis behaves similarly, P tsoo = P tsoo , b ( x ⊗ x + y ⊗ y ) + P tsoo , c z ⊗ z . The elastic dipole associated with the tetrahedral has two components, P t = P t , b ( x ⊗ x + y ⊗ y ) + P t , c z ⊗ z , corresponding to the basal P t , b and c -axis P t , c components. The tetrahedral-to-tetrahedraltransition along the c -axis also has hexagonal symmetry, P tstt = P tstt , b ( x ⊗ x + y ⊗ y ) + P tstt , c z ⊗ z . Finally, the octahedral-to-tetrahedral transition state does not locally have hexagonal sym-metry. If we consider one jump along a [1¯10 c a ], oriented so that it lies in the yz plane, thenwe can write the three components of the elastic dipole, P tsot = P tsot , b ( x ⊗ x + y ⊗ y ) + P tsot , t ( − x ⊗ x + y ⊗ y ) + P tsot , c z ⊗ z , corresponding to a basal, tetragonal, and c -axis contributions. Alternatively, we could writecontributions as P tsot , (cid:107) = ( P tsot , b + P tsot , t ) y ⊗ y along the hop and P tsot , ⊥ = ( P tsot , b − P tsot , t ) x ⊗ x perpendicular to the hop. We find that, after applying the 3-fold rotation axis that thetetragonal component only appears in our bare term.Hexagonal symmetry reduces the elastodi ff usion tensor to seven independent compo-nents, of which two have contributions from correlation. The contributions from correla-tion come from the tetrahedral site. The change in the velocity vector for the tetrahedralsite is δ b t = ± (cid:32) c λ t → o (cid:34) P tsot , b −
12 ( P t , b + (cid:104) P (cid:105) b ) (cid:35) − c λ t → t (cid:34) P tstt , b −
12 ( P t , b + (cid:104) P (cid:105) b ) (cid:35)(cid:33) ρ / ( z ⊗ x ⊗ x + z ⊗ y ⊗ y ) ± (cid:32) c λ t → o (cid:34) P tsot , c −
12 ( P t , c + (cid:104) P (cid:105) c ) (cid:35) − c λ t → t (cid:34) P tstt , c −
12 ( P t , c + (cid:104) P (cid:105) c ) (cid:35)(cid:33) ρ / z ⊗ z ⊗ z (56)15 uly 13, 2016 Philosophical Magazine paper and for the symmetric rate matrix δω = − (cid:16) λ t → t (cid:104) P tstt , b − P t , b (cid:105) + λ t → o (cid:104) P tstt , b − P t , b (cid:105)(cid:17) ( x ⊗ x + y ⊗ y ) − (cid:16) λ t → t (cid:104) P tstt , c − P t , c (cid:105) + λ t → o (cid:104) P tstt , c − P t , c (cid:105)(cid:17) z ⊗ z , (57)where (cid:104) P (cid:105) = P o ρ o + P t ρ t , as in our FCC lattice. Symmetry gives us 6 unique elastodi ff usioncomponents, as basal-plane isotropy requires 2 d = d − d . Moreover, only the d = d and d components have contributions from correlations; only d = d , d , d = d and d have contributions from the geometric terms; and d = d only has contributionsfrom geometry. Using the Voigt notation for fourth-rank tensors with Eqn. 34, d = d = a (cid:40) + β (cid:32) P tsot , b + P tsot , t − (cid:104) P (cid:105) b (cid:33)(cid:41) ρ o λ o → t d = c (cid:40) (cid:104) + β (cid:16) P tsoo , c − (cid:104) P (cid:105) c (cid:17)(cid:105) ρ o λ o → o + (cid:34) + λ t → o λ t → o + λ t → t (cid:16) P tstt , c − (cid:104) P (cid:105) c (cid:17) + λ t → t λ t → o + λ t → t (cid:16) P tsto , c − (cid:104) P (cid:105) c (cid:17)(cid:35) ρ t λ t → o λ t → t λ t → o + λ t → t (cid:41) d = d = a (cid:40) β (cid:32) P tsot , b − P tsot , t − (cid:104) P (cid:105) b (cid:33)(cid:41) ρ o λ o → t d = d = a (cid:110) β (cid:16) P tsot , c − (cid:104) P (cid:105) c (cid:17)(cid:111) ρ o λ o → t d = d = c (cid:40) β (cid:16) P tsoo , b − (cid:104) P (cid:105) b (cid:17) ρ o λ o → o + (cid:34) λ t → o λ t → o + λ t → t (cid:16) P tstt , b − (cid:104) P (cid:105) b (cid:17) + λ t → t λ t → o + λ t → t (cid:16) P tsto , b − (cid:104) P (cid:105) b (cid:17)(cid:35) ρ t λ t → o λ t → t λ t → o + λ t → t (cid:41) d = d = a ρ o λ o → t + c ρ o λ o → o + c ρ t λ t → o λ t → t λ t → o + λ t → t d = d − d = a (cid:26) + β P tsot , t (cid:27) ρ o λ o → t (58)and all other terms are zero. The activation volume tensor can be derived as before; thescalar activation volume from Eqn. 4 is V act = k B T B (cid:40) (cid:16) D − d + D − d + D − d (cid:17) S + S + S S + S + S + S + (cid:16) D − d + D − d (cid:17) S + S S + S + S + S (cid:41) , (59)where the bulk modulus B − = S + S + S + S . Carbon in iron
Carbon in body-centered cubic iron sits on octahedral sites, and transitions through tetra-hedral transition states. Veiga et al. [9] used an EAM potential to compute the octahe-dral and tetrahedral energies and elastic dipoles. The lattice constant for this potential16 uly 13, 2016 Philosophical Magazine paper is a = . C =
243 GPa, C =
145 GPa, and C =
116 GPa. The tetrahedral transition state is 0.816 eV above the octahedral site, andthe attempt frequency is taken as 10 THz[9]. The dipole tensors can be separated into par-allel and perpendicular components; the parallel direction points towards the closest Featoms for the C, while the perpendicular components lie in the interstitial plane. For theoctahedral, the parallel component is 8.03 eV, and the perpendicular components are both3.40 eV; for the tetrahedral transition state, the parallel component is 4.87 eV, and theperpendicular are both 6.66 eV. T − [K − ] -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 D [ c m / s ] Dd d d = D T [K] Figure 1. Di ff usivity and elastodi ff usivity of C in BCC Fe between 300K and 1200K. The carbon di ff usivity is Arrhenius,with D ( T ) = . × − exp( − . / k B T ) cm / s. The derivative d = dD xx / d ε xx changes sign: negative below ≈ d = dD xy / d ε xy = D , as it only has a contribution from the geometric change with strain. T − [K − ] V [ a t o m i c v o l u m e ] V total V V V T [K] Figure 2. Activation volume of C in BCC Fe. The total activation volume from Eqn. 4 is the sensitivity to the isotropicdi ff usivity to hydrostatic pressure; at 300K it is 0 . Ω Fe = .
357 Å . The components from Eqn. 3 have cubic symmetry,and so the representative terms at 300K are V = V xxxx = − . Ω Fe = − .
368 Å , V = V xxyy = . Ω Fe = .
220 Å ,and V = V xyxy = . Ω Fe = .
142 Å . The only contribution to V xyxy is geometric. Fig. 1 shows the Arrhenius behavior of the di ff usivity, while Fig. 2 reveals unusual be-havior in the activation volume from 300K to 1200K. The di ff usivity follows a simpleArrhenius relationship. The elastodi ff usion tensor has three unique components: d = dD xx / d ε xx , d = dD xx / d ε yy , d = dD xy / d ε xy . The d term only contains geometriccontributions, while d changes sign as the (positive) geometric contribution dominates athigh temperatures and the (negative) activation contribution dominates at low temperatures.The d term is an order of magnitude larger than the other elastodi ff usivity coe ffi cients,17 uly 13, 2016 Philosophical Magazine paper and also leads to the unusual activation volume behavior where V is negative, V is pos-itive, and V total = V + V is positive. This means that for a compressive uniaxial stressalong (cid:104) (cid:105) , the di ff usivity for carbon increases along (cid:104) (cid:105) , while the perpendicular dif-fusivities decrease . The unusual behavior is caused by the Poisson e ff ect combined withthe large d value and small d value. This still occurs while a hydrostatic compressivestress decreases the di ff usivity.
4. Conclusions
The general expression for the elastodi ff usion tensor, activation energy and volume tensorsbuild on previous work that assumed interstitial migration without correlation. As shownhere, correlation contributions must be included even in cases of interstitial di ff usion wheredictated by the di ff usion network. The application of this method extends beyond intersti-tial solutes in a crystal, but can also be used for more complex cases such as self-interstitials(which may possess multiple intermediate states) or di ff usion through boundaries. There isan open-source numerical implementation of the di ff usion, elastodi ff usion, and activationbarrier tensors, described in Appendix C[19]. The perturbation analysis laid out here canalso be applied to error propagation, by calculating the derivative of di ff usivity with respectto site and transition energies. It should be possible, though more complicated, to derivesimilar expressions for the elastodi ff usion tensor for more complex di ff usion mechanisms,such as vacancy-mediated solute di ff usion. Acknowledgements
The author thanks Maylise Nastar and Pascal Bellon for helpful conversations.
Disclosure
The author has no conflicts of interest.
Funding
Research supported in part by the U.S. Department of Energy, O ffi ce of Basic En-ergy Sciences, Division of Materials Sciences and Engineering under Award ffi ce of Naval Research grant N000141210752, and in part by the National Science Foun-dation Award 1411106. Part of this research was performed while the author was visitingthe Institute for Pure and Applied Mathematics (IPAM) at UCLA, which is supported bythe National Science Foundation (NSF). References [1] C.P. Flynn,
Point defects and di ff usion , Clarendon Press, Oxford, 1972.[2] A.R. Allnatt and A.B. Lidiard, Atomic transport in solids , chap. 5, Cambridge UniversityPress, Cambridge (1993), p. 202–203. uly 13, 2016 Philosophical Magazine paper [3] H.H. Wu and D.R. Trinkle, Direct di ff usion through interpenetrating networks: Oxygen intitanium , Phys. Rev. Lett. 107 (2011) 045504.[4] F.C. Larche and J.W. Cahn, The e ff ect of self-stress on di ff usion in solids , Acta metall. 30(1982) 1835–1845.[5] M.J. Aziz, Thermodynamics of di ff usion under pressure and stress: Relation to point defectmechanisms , Appl. Phys. Lett. 70 (1997) 2810–2812.[6] M.S. Daw, W. Windl, N.N. Carlson, M. Laudon and M.P. Masquelier, E ff ect of stress on dopantand defect di ff usion in Si: A general treatment , Phys. Rev. B 64 (2001) 045205, Available at http://link.aps.org/doi/10.1103/PhysRevB.64.045205 .[7] M.J. Aziz, Y. Zhao, H.J. Gossmann, S. Mitha, S.P. Smith and D. Schiferl, Pressure and stresse ff ects on the di ff usion of B and Sb in Si and Si-Ge alloys , Phys. Rev. B 73 (2006) 054101.[8] R.G.A. Veiga, M. Perez, C.S. Becquart, C. Domain and S. Garruchet, E ff ect of the stress fieldof an edge dislocation on carbon di ff usion in α -iron: Coupling molecular statics and atomistickinetic Monte Carlo , Phys. Rev. B 82 (2010) 054103.[9] R.G.A. Veiga, M. Perez, C. Becquart, E. Clouet and C. Domain, Comparison of atomistic andelaticity approaches for carbon di ff usion near line defects in α -iron , Acta mater. 59 (2011)6963–6974.[10] C.H. Woo and C.B. So, The e ff ect of stress on point-defect di ff usion in hcp metals and irradi-ation creep , Philos. Mag. A 80 (2000) 1299–1318.[11] W.L. Chan, R.S. Averback and Y. Ashkenazy, Anisotropic di ff usion of point defects in metalsunder a biaxial stress field simulation and theory , J. Appl. Phys. 104 (2008) 023502.[12] E.J. Savino, Point defect-dislocation interaction in a crystal under tension , Philos. Mag. 36(1977) 323–340.[13] P.H. Dederichs and K. Schroeder,
Anisotropic di ff usion in stress fields , Phys. Rev. B 17 (1978)2524–2536, Available at http://link.aps.org/doi/10.1103/PhysRevB.17.2524 .[14] H. J´onsson, G. Mills and K.W. Jacobsen, Nudged elastic band method for finding minimumenergy paths of transitions , in
Classical and Quantum Dynamics in Condensed Phase Sim-ulations , B.J. Berne, G. Ciccotti and D.F. Coker, eds., World Scientific, Singapore, 1998, p.385–404.[15] G. Henkelman, B.P. Uberuaga and H. J´onsson,
A climbing image nudged elastic band methodfor finding saddle points and minimum energy paths , J. Chem. Phys. 113 (2000) 9901–9904.[16] A. Janotti, M. Krˇcmar, C.L. Fu and R.C. Reed,
Solute di ff usion in metals: Larger atoms canmove faster , Phys. Rev. Lett. 92 (2004) 085901.[17] M. Mantina, Y. Wang, L.Q. Chen, Z.K. Liu and C. Wolverton, First principles impurity di ff u-sion coe ffi cients , Acta mater. 57 (2009) 4102–4108.[18] G.H. Vineyard, Frequency factors and isotope e ff ects in solid state rate processes , J. Phys.Chem. Solids 3 (1957) 121–127.[19] D.R. Trinkle, O nsager , http://dallastrinkle.github.io/Onsager (2016),doi: // / zenodo.57407.[20] R. Agarwal and D.R. Trinkle, Light element di ff usion in Mg using first principles calculations:Anisotropy and elastodi ff usion , 2016. arXiv:1605.05650.[21] M. Glazer and G. Burns, Space Groups for Solid State Scientists , 3rd ed., Elsevier, 2013.[22] Y. Hanlumyuang, P. Gordon, T. Neeraj and D. Chrzan,
Interactions between carbon solutesand dislocations in bcc iron , Acta mater. 58 (2010) 5481–5490, Available at .[23] C. Varvenne, F. Bruneval, M.C. Marinica and E. Clouet,
Point defect modeling in materials:Coupling ab initio and elasticity approaches , Phys. Rev. B 88 (2013) 134102, Available at http://link.aps.org/doi/10.1103/PhysRevB.88.134102 .[24] T. Garnier, V.R. Manga, P. Bellon and D.R. Trinkle, Di ff usion of Si impurities in Ni understress: A first-principles study , Phys. Rev. B 90 (2014) 024306.[25] H. Kim and D.R. Trinkle, Size misfits of point defects in aliovalently doped SrTiO usingdensity functional theory , Comp. Mat. Sci. 119 (2016) 41–45.[26] T. Hahn (ed.), International Tables for Crystallography , 4th ed., Vol. A, Dordrecht: KluwerAcademic, 1996. uly 13, 2016 Philosophical Magazine paper [27] S.C. Middleburgh and R.W. Grimes, Defects and transport processes in beryllium , Acta mater.59 (2011) 7095–7103.[28] P. Zhang, J. Zhao and B. Wen,
Retention and di ff usion of H, He, O, C impurities in Be , J. Nucl.Mater. 423 (2012) 164–169.[29] A. O’Hara and A.A. Demkov, Oxygen and nitrogen di ff usion in α -hafnium from first prin-ciples , Appl. Phys. Lett. 104 (2014), 211909, Available at http://scitation.aip.org/content/aip/journal/apl/104/21/10.1063/1.4880657 .[30] G.W. Stewart, On the perturbation of pseudo-inverses, projections and linear least squaresproblems , SIAM Review 19 (1977) 634–662, Available at .[31] L.W. Johnson, R.D. Riess and J.T. Arnold,
Introduction to Linear Algebra , chap. 3, Pearson(2002), p. 199–201.
Appendix A. Pseudoinverse perturbation
The expression in Eqn. 21, which is similarly used to express the first-order change inpseudoinverse of ω with respect to strain, comes from the perturbation of the pseudoinverseand lacks terms that would correspond to the perturbative change in the null space. In thiscase, the simplified expression is due to the symmetry of the matrix ω , the perturbativechanges δω , and the symmetry of the vector dot-products with b . If a general matrix A haspseudoinverse A + , then we can expand the pseudoinverse B + = ( A + E ) + up to first orderin the perturbation E as B + = A + − A + ( AA + ) E ( A + A ) A + + ( A † A ) + ( A + A ) E † ( − AA + ) − ( − A + A ) E † ( AA + )( AA † ) + + O ( (cid:107) E (cid:107) )(A1)from Eqn. (3.24) of [30]. For our case of interest, A = ω and A + = g are both real symmetricmatrices; moreover, we want to compute v · B + · v for an arbitrary vector v , v · B + · v = v · A + · v − v · ( A + EA + ) · v + v · (cid:26) ( A ) + ( A + A ) E ( − AA + ) − (cid:104) ( A ) + ( A + A ) E ( − AA + ) (cid:105) T (cid:27) · v + O ( (cid:107) v · E } ) = v · A + · v − v · ( A + EA + ) · v + O ( (cid:107) v · E } ) (A2)which, for our problem, reduces to first order in δω b · δ g · b = − b · g δω g · b = − γ · δω · γ . (A3) Appendix B. Vector and tensor bases at sites
To construct the vector and tensor bases, we first must determine the eigenvectors of eachgroup operation. There are only ten types of operations available in a crystal: identity, a 2-,3-, 4-, or 6-fold axis, mirror through a plane, or a 2-, 3-, 4-, or 6-fold axis combined witha mirror through that same plane. These can be identified from the determinant—which is20 uly 13, 2016 Philosophical Magazine paper − − − + mirror1 1 4-fold axis − − + mirror1 0 3-fold axis − − + mirror1 − − − e , e , e —inversion—triply-degenerate eigenvalues of −
1, for the three eigenvectors e , e , e —andmirror—one eigenvalue of − v corresponding to the mirror plane normaland double-generate eigenvalues of 1 for the two orthonormal vectors v (cid:48) and v (cid:48)(cid:48) that areorthogonal to v . For the remaining cases, the rotational axis v is always an eigenvector withan eigenvalue + − n th roots of unity e ± i π/ n for the n -fold rotation axis. If v (cid:48) and v (cid:48)(cid:48) are two orthonormal vectors that are also orthogonal to the rotational axis v , thenthe eigenvectors are ( v (cid:48) ± i v (cid:48)(cid:48) ) / √ d -ranktensors can be constructed from di ff erent dyad products of eigenvectors of the group oper-ation. In particular, the product of the eigenvalues of the eigenvectors needs to be 1. Hence,for a vector basis (a 1-rank tensor), we have: (a) identity produces a 3-dimensional basisof { e , e , e } , (b) mirror produces a 2-dimensional basis of { v (cid:48) , v (cid:48)(cid:48) } , two orthonormal vec-tors in the mirror plane, (c) n -fold rotational axes without mirror produces a 1-dimensionalbasis of { v } , the rotational axis vector, and (d) all other cases produce a 0-dimensional nullbasis, ∅ . Each group operation produces its own basis, and these vector bases are then in-tersected to produce the final basis, using some simple rules: if p - and q -dimensional basesare intersected, the resulting dimension of the final basis is not greater than min { p , q } . Then,the remaining non-trivial cases are: two 2-dimensional, a 2- and 1-dimensional, and two1-dimensional intersections. In the first case, if v and v (cid:48) are the vectors made by the cross-product of the two vectors in each basis, then the resulting intersection will be either (a)the same 2-dimensional basis if v is parallel to v (cid:48) , or (b) a 1-dimensional basis consistingof v × v (cid:48) , otherwise. In the second case, if v is the basis for the 1-dimensional basis, and v is the cross-product of the two vectors in the 2-dimensional basis, then the resultingintersection will be either (a) 1-dimensional basis consisting of v if v is perpendicularto v , or (b) the 0-dimensional empty basis ∅ , otherwise. Finally, if v and v (cid:48) are the two1-dimensional bases, the resulting intersection will either be (a) v if v is parallel to v (cid:48) , or(b) the 0-dimensional empty basis ∅ , otherwise.The construction of the second-rank tensor bases is more complicated, but reduces tothree distinct cases. We are specifically interested in symmetric second-rank tensors, so thelargest the bases can be is 6-dimensional. For identity and inversion, the full 6-dimensionalbasis is invariant: { e ⊗ e , e ⊗ e , e ⊗ e , ( e ⊗ e + e ⊗ e ) / √ , ( e ⊗ e + e ⊗ e ) / √ , ( e ⊗ e + e ⊗ e ) / √ } . For a mirror operation or a 2-fold rotation, if v is the mirror plane normal / rotational axis, and v (cid:48) and v (cid:48)(cid:48) are the remaining orthonormal vectors, the 4-dimensionalbasis is invariant: { v ⊗ v , v (cid:48) ⊗ v (cid:48) , v (cid:48)(cid:48) ⊗ v (cid:48)(cid:48) , ( v (cid:48) ⊗ v (cid:48)(cid:48) + v (cid:48)(cid:48) ⊗ v (cid:48) ) / √ } . Finally, for the n -fold rota-tions (with or without a mirror), if v is the rotational axis, and v (cid:48) and v (cid:48)(cid:48) are the remaining21 uly 13, 2016 Philosophical Magazine paper orthonormal vectors, the 2-dimensional basis is invariant: { v ⊗ v , ( v (cid:48) ⊗ v (cid:48) + v (cid:48)(cid:48) ⊗ v (cid:48)(cid:48) ) / √ } .To intersect any two bases, the simplest approach is to find the null-space of the combinedcolumn spaces constructed from the bases using singular-value decomposition[31]. Appendix C. Implementation
A full numerical implementation of the algorithms in Python described are available ongithub[19] under the MIT License. This includes algorithms to analyze a given crystal(lattice and atomic basis), find generators for the space group operations, determine allpoint group operations for each site, identify Wycko ff positions, generate vector and ten-sor bases, construct a jump network in the crystal, and identify unique jumps. Once theenergies, prefactors, and elastic dipoles are determined for the unique sites and jumps, thenumerical implementation can compute the di ff usivity, elastodi ff usion, and activation bar-rier tensors for a given temperature. The implementation includes numeric examples ofthe analytic results from this paper, with input files for FCC, BCC, and HCP interstitialdi ffff