Non-Gaussian diffusion profiles caused by mobile impurity-vacancy pairs in the five frequency model of diffusion
aa r X i v : . [ c ond - m a t . m t r l - s c i ] J a n Non-Gaussian diffusion profiles caused by mobile impurity-vacancy pairs in the fivefrequency model of diffusion
V. I. Tokar IPCMS, Universit´e de Strasbourg–CNRS, UMR 7504, 23 rue du Loess, F-67034 Strasbourg, France (Dated: October 12, 2018)Vacancy-mediated diffusion of impurities under strong impurity-vacancy (I-v) attraction has beenstudied in the framework of the five-frequency model (5FM) for the FCC host. The system ofimpurities and tightly bound I-v pairs has been treated in the framework of the rate-equationsapproach of Cowern et al., Phys. Rev. Lett. , 2434 (1990), developed for the description of thenon-Gaussian diffusion profiles (NGDPs) observed in dopant diffusion in silicon. In the present studythis approach has been extended to derive a three-dimensional (3D) integro-differential equationdescribing the pair-mediated impurity diffusion. The equation predicts the same 1D NGDPs as inCowern et al. but can be also used for the simulation of 3D profiles of arbitrary geometry in thesystems where the diffusion proceeds via a mobile state. The parameters of the theory has beencalculated within the 5FM on the basis of available literature data. The database on impuritiesin aluminum host has been analyzed and promising impurity-host systems for the observation ofNGDPs has been identified. The diffusion profiles for an impurity where NGDPs are expected to beeasily detectable have been simulated. It has been argued that with the input parameters calculatedon the basis of experimental diffusion constants the simulated NGDPs can be accurate enough toserve as a quantitative test of the 5FM. I. INTRODUCTION
Miniaturization of electronic devices to nanometersizes has necessitated investigation of peculiarities thattechnologically important processes may exhibit at thisscale.
One such process that plays a key role in the dop-ing of semiconductor chips is the diffusion of impurities inelemental crystalline hosts. In experiments of Cowern etal. it was established that some dopants in the siliconhost exhibit exponential diffusion profiles instead of theconventional Gaussian ones. This seemingly non-Fickianbehavior was attributed to the diffusion mediated by amobile intermediate state of impurities induced by inter-action with point defects, such as vacancies and inter-stitials. A phenomenological model based on the notionof the mobile state developed by Cowern et al. sat-isfactorily described the experimental observations withthe use of only two adjustable parameters. However,the development of the model into a quantitative the-ory has been hampered by insufficient understanding ofthe microscopic mechanisms underlying the mobile state.The main problem poses the presence in semiconductorhosts of several competing mechanisms with their rel-ative importance being impurity-specific and not fullyunderstood.
The analysis of defect-assisted diffusion considerablysimplifies when only one mechanism dominates, as wasthe case in the diffusion in the Cu(001) surface layerwhere exponential tails in diffusion profiles were alsoobserved.
The vacancy mechanism was shown to bedominant and due to its simplicity it was possiblenot only to reliably fit the values of microscopic param-eters to experimental data but also to confirm them inthe first-principles calculations. Diffusion in two dimen-sions (2D), however, is qualitatively different from 3Dcase and so cannot serve as a model of vacancy- mediated diffusion in 3D systems.But the vacancy mechanism is also common in 3Dsolids, even in semiconductors.
Arguably, it has beenbest studied in the FCC metals in the framework of theclassic five frequency model (5FM).
There exists awealth of literature on the pertinent host-impurity sys-tems, on approximate solutions of the model, on compari-son with experiment as well as large databases containingsystematized data on both first-principles calculations ofthe model parameters and on their empirical values (see,e. g., Refs. 21–32 and references therein).Furthermore, the mobile state of impurity in the 5FMis also known.
It appears in the case of strongimpurity-vacancy (I-v) attraction when the I-v pairs thatform can exercise long sequences of diffusion jumps be-cause in the bound state the diffusion-mediating vacancyis permanently available. Similarly to the 2D case dis-cussed above, this mechanism was shown to pro-duce non-Gaussian diffusion profiles (NGDPs) also in theFCC systems. The term NGDPs will be used through-out the paper to refer to the profiles that Cowern et al.called exponential because unlike in the surface layerwhere truly exponential behavior can be observed due tothe specifics of the STM experiment, below we willsee that in 3D the tails of the profiles are asymptoticallyalways Gaussian with the exponential behavior being ob-servable only at intermediate distances.The aim of the present paper is to develop a quan-titative theory of impurity diffusion propagated by thebound I-v pairs using the framework of the 5FM for FCChosts. The theory will be essentially based on the phe-nomenological approaches of Refs. 4 and 13 with the in-put parameters calculated within the 5FM from the dataavailable in literature sources.
The paper is organized as follows. In the next sectionthe 5FM is briefly introduced; in Sec. III theory of dif-fusion of individual tightly-bound I-v pairs is presented;in Sec. IV the diffusion of the ensemble of immobile im-purities and mobile pairs is treated within a phenomeno-logical theory based on the approach of Cowern et al. ;in Sec. V the integro-differential equation describing 3Ddiffusion of impurities mediated by the mobile state isderived and illustrated by simulation of 1D NGDPs forlanthanum impurity in aluminum host; in the concludingsection VI the results obtained are briefly summarized. II. THE 5FM
The 5FM for the FCC lattice is a representativeof the class of models that describe the vacancy-mediateddiffusion as a stochastic process characterized by a set ofthe transition rates or frequencies of the vacancy jumpsbetween the sites of the host lattice. In Fig. 1 the mean-ing of the five frequencies w k , k = 0 −
4, for the FCC hostare explained. In the canonical model that will be usedin the present paper the vacancy jumps are restricted toonly nearest neighbor (NN) sites though generalizationson more complex models with larger sets of parametersare possible.
As can be seen from the defini-tion of 5FM in Fig. 1, in the vicinity of impurity onlyone vacancy is assumed to be present. This, of course,is an approximation but it will be sufficient for our pur-poses because the NGDPs we are going to study are themost pronounced at low temperatures where the vacancyconcentration is very small. For example, using the ex-perimental vacancy formation enthalpy in aluminum E f = 0 .
67, the vacancy concentration at temperature50 ◦ C that will be used in our simulations in Sec. V canbe estimated to be c v ≃ e − E f /k B T = 3 . · − . (1)Here we neglected the entropic contribution (in a rigoroustreatment the Gibbs free energy should be used instead ofthe energy) because the formation entropies s f are usu-ally quite small s f ≈ and the temperature 50 ◦ C inenergy units is ∼ .
03 eV while the errors in both experi-mental and theoretical definitions of E f are considerablylarger being at least of order of 0.1 eV. The jump frequencies can be calculated from the valuesof the activation barriers E k and the attempt frequencies ν k (see Table I) as w k = ν k e − E k /k B T . (2)Because the stochastic dynamics in the 5FM is governedby thermal excitations, it satisfies the detailed balancecondition which establishes the following relation be-tween the frequencies and the binding energy E b of theI-v pair w w ≃ exp (cid:18) − E b k B T (cid:19) . (3) w w w w FIG. 1. Jump frequencies of the vacancy (gray square) inthe vicinity of the impurity (black circle) as defined in the5FM: w –the jumps in the first coordination shell (CS) of theimpurity; w –frequency of exchange with the impurity atom; w –frequency of dissociative jumps away from the impurityinto higher CS; w –associative jumps from higher CS intothe first one. Not shown in the figure is the frequency w ofexchange with atoms in the host bulk.TABLE I. Parameters entering Eqs. (2) and (9) correspondingto lanthanum impurity (I) in aluminum host (H) taken fromthe database of Ref. 31 k E k (eV) ν k (THz)0 0.5814 4.22421 1.18 3.98682 0.0623 2.51823 0.5637 4.19754 0.092 4.0664 Q (eV) D (cm /s)I 0.7776 0.081112H 1.2661 0.064623 In this definition E b is assumed to be positive for I-vattraction. As in Eq. (1), here we also neglected the en-tropy contribution because of its smallness in comparisonwith the errors in E b . The frequencies w – w together with the vacancy con-centration from Eq. (1) and the value of the lattice pa-rameter a can be used to calculate the diffusion constantof the impurity in the 5FM as D = c v w w w f a , (4)where the correlation factor f = 2 w + 7 F ( w /w ) w w + 2 w + 7 F ( w /w ) w . (5)Accurate expressions for F were derived in Refs. 41 and42. For our purposes will be necessary the values of f inthe case of strong I-v binding (large E b ). According toEq. (3) this corresponds to cases when w → w → ∞ . From Eq. (5) it is seen that in the first case( w →
0) the value of F is irrelevant while in the secondcase it is known exactly: F ( w /w → ∞ ) = 2 / . (6)Thus, the correlation factor in the limit of infinitelystrong binding is f ∞ = w + w w + w + w . (7)Eq. (5) also accurately reproduces the exactly knownvalue for the case of self-diffusion when all w i are equal: f ≃ . . (8)The most serious obstacle to quantitative predictionsbased on the 5FM is that the parameters of the modeleither as fitted to experimental data or as obtained infirst-principles calculations may contain quite significanterrors. The problem aggravates at lower tem-peratures because of the Arrhenius behavior in Eq. (2).In order to reconcile the measured and calculated valuesof th host ( H ) and the impurity ( I , though this super-script will be usually omitted for brevity) diffusion con-stants D H,I = D H,I exp (cid:18) − Q H,I k B T (cid:19) , (9)a correction coefficient was suggested in Ref. 30 C shift = A shift exp (cid:18) − E shift k B T (cid:19) . (10)In the case of aluminum A shift = 12 and E shift = 0 . ◦ C the coefficient C shift ≈ − which meansthat the discrepancy between the theory and the exper-iment amounts to two orders of magnitude at this tem-perature. Because various quantities in the 5FM maydepend on six parameters (five frequencies plus the va-cancy formation energy) and their errors combine, wewill try to maximally reduce their number in our calcu-lations and whenever possible use instead experimentallymeasurable quantities. III. DIFFUSION OF TIGHTLY BOUND I-VPAIRS
According to Eq. (3), strong I-v binding meaning large E b value takes place when either w is very small or w is very large or both. The case of vanishing escape fre-quency w was solved in Ref. 19 so we will consider amore general case of large w ≫ w which comprises alsothe small w case.To begin with, let us find the solution in the limitof the infinitely large of w and for definiteness let usrestrict our consideration only to diffusion along Z di-rection because in cubic lattices all h i directions areequivalent. Similar to Ref. 19, the problem in this casecan be reduced to the solution of a set of three equationsfor the probabilities of three inequivalent mutual I-v ori-entations. They correspond to the vacancy being in oneof three classes 1–3 of the sites in the first coordination ✻ Za − a ✧✧✧✧✧✧✧✧✧✧✧✧ ✧✧✧✧✧✧✧✧✧✧✧✧ ✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧ ⑦✉ ✉ ✉ ✉ ✉ ✉ ✉ ✉ ✉ ✉ ✉ ✉ ✉ ✉ ✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧ ✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧✧ +3 +3+1 +1 ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ ❡ FIG. 2. Axisymmetric diffusion along Z direction of the FCClattice as described by the 5FM: Distribution of lattice sitesin the vicinity of the impurity (black circle) among 13 equiv-alence classes. Black points—the classified sites in the cubesvertexes; crosses—the face centered sites on the visible facesin the drawing, open circles—on the invisible faces. shell (CS) of the impurity (see Fig. 2); inside these classesall positions are equivalent due to our choice of the sym-metry direction. In Fig. 2 all sites pertinent to the 5FMare divided into 13 equivalence classes. In the case ofinfinitely large w the vacancy will spend all its time onthe sites in classes 1–3 because of the following. Classesfrom 4 to 13 comprise the sites that can be reached fromthe first CS in one jump with the rate w . If the vacancyjumps at one of these sites, it will have at least one NNsite in the first CS (the one it just jumped from) but ingeneral there will be l ≥ lw or diffuses further away from the impurity withthe rate (12 − l ) w . The probability of jumping back tothe first CS is given by the ratio p back = lw lw + (12 − l ) w (11)which tends to unity when w → ∞ . Moreover, becausethe residence time ∆ t of the vacancy at the “outer” siteis inversely proportional to the total rate in the denomi-nator of Eq. (11) ∆ t ∝ / [ lw + (12 − l ) w ], as w → ∞ the time goes to zero. Thus, the vacancy spends all itstime in the first CS of impurity with the excursions tosites in classes 4-13 serving only to jumps between classes1-3 in the first CS.To formalize this picture let us number the (001) planesin the direction of the diffusion by integer numbers m, n ranging from minus to plus infinity. When the impurityis positioned in plane n the sites of class 1 are in theplane m = n + 1, sites of class 2 are in the same plane( n ) and those of class 3 in the plane m = n − W αα ′ nm as ddt C αn ( t ) = X m,α ′ ˜ W αα ′ nm C α ′ m ( t ) , (12) where C αn is the concentration of pairs on plane n withthe vacancy being in class α = 1 −
3. The transitionmatrix that has the standard loss-gain structure is˜ W nm = − (2 w + w + w ) δ nm (2 w + w ) δ nm w δ m − n, + w δ nm (2 w + w ) δ nm − (4 w + 3 w ) δ nm w + w ) δ nm w δ n − m, + w δ nm (2 w + w ) δ nm − (2 w + w + w ) δ nm , (13)where the rows and columns are numbered by the classindexes α, α ′ while subscripts m and n denote the planes.The non-diagonal entries of the matrix correspond totransition rates between the classes. For example, ˜ W nm is the rate of the vacancy jumps from a site in class 1 to asite in class 2 as well as a possible ensuing displacementof the impurity from plane m to plane n . The jump canproceed by several routes. The vacancy can reach class 2directly via NN jumps. This contributes 2 w to the ratebecause there are exactly two sites in class 2 that are NNto a site from class 1. The indirect way is to first jumpout of the first CS on one of two NN sites in class 8 oron the site in class 5 (in total three possibilities). Butthe return jumps will end up in class 2 only in half of thecases which reduces the contribution to 3 w /
2. The ma-trix element is proportional to δ mn because the vacancymoves do not change the impurity position. In fact, onlydirect I-v exchanges lead to the impurity diffusion, so theonly contributions that change the position of the impu-rity are those proportional to w , which is reflected inEq. (13) (see Fig. 1).To simplify the task of solving Eq. (12), let us restrictour attention only to macroscopic diffusion that developsat a spatial scale much larger than the lattice constant a .At this scale the positions of the impurity and the nearbyvacancy are essentially the same so concentrations C αn need not be distinguished and only evolution of the totalpair density C n = X α =1 C αn (14)will be detectable experimentally via the impurity pro-files. Long diffusion distance means long duration. Butsolutions of stochastic evolution equations of the kind ofEq. (12) tend toward the equilibrium via exponentiallyattenuating modes ∝ exp( z i t ) with all z pi being less orequal to zero. So we will be interested only in the longest-lived modes that correspond to the smallest z i .For a translationally-invariant system with constantcoefficients the standard way to find the attenuation ratesis to reduce differential equation Eq. (12) to an algebraicequation with the use of integral transforms which wechose to be the Laplace transform in time variable andthe Fourier transform in the spatial variables (the LFtransform). The 1D Fourier transform of the transition matrix Eq. (13) is˜ W K = − x − y x y + γ K x − x xy + γ ∗ K x − x − y , (15)where x = 2 w + 32 w (16) y = w + 14 w (17) γ K = w ( e iaK/ −
1) (18)and a/ t with the pa-rameter z we obtain the characteristic equationdet( z − ˜ W K ) = 0 . (19)As explained above, we do not need exact expressions forthe solutions of this equation but are interested only inthe smallest eigenvalue that describes the pair diffusion.At small Fourier momenta K corresponding to large dis-tances in real space the eigenvalue will tend to zero as ∼ K because the matrix in Eq. (15) is Hermitian. Thesmallest root (let us denote it z ) can be easily found tothis accuracy from Eq. (19) from the linearized equationbecause as z → z can be dropped: z p ( K ) ≃ w ( w + w )12( w + w + w ) a K . (20)Here superscript p reminds us that the solution describesthe pair diffusion and the momentum (0 , , K ) describingthe diffusion along Z direction is replaced by a generalmomentum K = ( K X , K Y , K Z ) because in cubic crystalsthe diffusion is isotropic to this order in K . The remain-ing two eigenvalues are finite at small K and up to terms O ( K ) are z p ( K ) = − x + O ( K ) ≃ − w − . w z p ( K ) = − x − y + O ( K ) ≃ − w + w + w )(21)as can be easily verified by direct substitution in Eq. (19)at K = 0. As is seen, the eigenvalues 2–3 remain finiteas K →
0. Therefore, only the term corresponding to z p ( K ) = O ( K ) will survive at large times so the Fouriertransform of the pair density Eq. (14) will behave as C K ( t ) | t →∞ ∝ exp [ − z p ( K ) t ] . (22)Differentiating this by t and taking the inverse Fouriertransform with respect to the spatial variables bringabout the conventional diffusion equation ∂C ( R , t ) /∂t ≃ D m ∇ C ( R , t ) , (23)where C is the continuum approximation to C n and ac-cording to Eqs. (20), (22), and (7) the diffusion constantof the I-v pairs is D m = w ( w + w ) a w + w + w ) = w f ∞ a . (24)For w = 0 this expression coincides with the result ofRef. 19.The most important property of D m is that in con-trast to the impurity diffusion constant Eq. (4), it is notproportional to the vacancy concentration and so can bemuch larger than D . Assuming strong binding in whichcase one can use f ∞ as f in Eq. (4) with the use of Eqs.(24), (1), and (3) one gets D m D ≃ w c v w ≃ e ( E f − E b ) /k B T . (25)Simple bond-counting arguments suggest that the va-cancy formation energy E f should be larger than thebinding energy E b . The latter can be found by compar-ing the energies of a vacancy surrounded by only hostatoms and the vacancy with one neighbor replaced bythe impurity which amounts to the difference betweenenergies of a single atomic bond. Creation of a vacancy,on the other hand, costs about six atomic bonds in theFCC lattice. So as T → ◦ C may reach according to Eqs. (1)–(3) and Table Ithe value ( D m /D ) T =50 ◦ C ≈ . (26)It is the presence of two modes of impurity diffusion withvery different diffusion constants that underlies the phe-nomenon of NGDPs. A. Decay rate of the I-v pair
Despite large diffusivity, the bound I-v pairs cannotdiffuse too far from the place of their association becauseof their finite lifetime. Irrespective of how strong the I-vbinding is, from Eq. (2) it is seen that at finite tempera-ture neither w can be strictly equal to zero nor w can beinfinitely large. Thus, if the lifetime of the pair is equal to τ decay its decay rate r = 1 /τ decay and the characteristicdistance λ of the pair diffusion before its decay is equal(up to a numerical constant) to the diffusion length ofthe pair during its lifetime λ = p D m τ decay = p D m /r. (27)In Refs. 4 and 5 this quantity is called the mean projectedpath length.In the 5FM the decay rate r can be found as the rateof definite separation of the vacancy from the impurity.It is often approximated by the rate of escape from thefirst CS into higher coordination shells which is equal to7 w because there is seven NN sites to a site in thefirst CS that do not belong to the first CS, as can beseen from Figs. 1 or 2. However, as we already saw, thisapproximation can be completely misleading in the caseof large w /w ratio when the vacancy returns into thefirst CS with probability that almost equals to unity. Toaccount for this r can be represented as a product of the“bare” escape rate 7 w and the renormalization factor p ∞ equal to the probability of definite I-v dissociationwhen the vacancy diffuses infinitely far away from theimpurity r = 7 w p ∞ . (28)The problem of calculating the vacancy escape prob-ability p ∞ is equivalent to finding the return probabil-ity into the first CS which is a standard problem of therandom walk theory (see, e. g., Ref. 46 and referencestherein). But to avoid complicated combinatorial cal-culations we will assess p ∞ with the help of numericalsimulations. The escape of the vacancy from the first CSinto the space beyond sufficiently large radius R max wassimulated for several values of the ratio w /w and fortwo radii R max = 100( a/
2) and 300( a/
2) and then inter-polated to R max = ∞ . The results are presented in Fig.3 together with the exactly known two end point values p ∞ ( w /w = 0) = 1 and p ∞ ( w /w = ∞ ) = 0 and withthe approximating interpolating expression p ∞ ≃ (1 + bw /w ) − , (29)where b ≈ .
35 was found from the largest simulatedratio w /w = 10 because large w /w values are themost interesting to our purposes. However, because ofthe diminishing number of the decays as w /w → ∞ ,the statistics are difficult to gather when the ratio is inthe range O (10 − ) that we are interested in (see thenext paragraph). It can be calculated in this case withinthe rigorous approach of Ref. 47, as will be shown in thesubsequent publication. Therefore, in our calculationsbelow we will use this more accurate value b = 1 . λ of the theory of Ref. 4 in terms of thefrequencies of the 5FM. To calculate λ for 50 impuritieslisted in the database of Ref. 31 for aluminum host use p ∞ w /w → ← w /w FIG. 3. Probability for the vacancy to definitely leave the im-purity neighborhood after escape from the first CS: bullets—the KMC simulation data and the exactly known end points;solid line: Eq. (29).TABLE II. Impurities in aluminum host with λ > a attemperature T = 50 ◦ C calculated from Eqs. (32) and (9)with the parameters taken from Ref. 31; E b = E − E .I λ (nm) E b (eV)La 328 0.47Se 161 0.44Te 34 0.46Bi 30 0.41Nd 22 0.27Ce 21 0.26Pb 17 0.36Tl 8.2 0.31Sb 4.5 0.30 has been made of Eqs. (27), (24), (28) and (29). Nine per-spective impurities with the largest values of λ were iden-tified (see Table II). It turned out that in all nine casesthe ratio w /w were very large, of order O (10 − ).The impurities with large λ are the most interestingones from experimental perspective, so it is desirable tocalculate the parameter with maximum precision. But,as was discussed at the end of the previous section, theaccuracy of the 5FM frequencies values is currently rathermodest and the errors may compound in the uncertaintyof the value of λ . So it would be preferable to express thelatter in terms of more reliable parameters. In the case w → ∞ this can be achieved as follows. At large w p ∞ | w ≫ w ≃ w bw . (30)Substituting this into Eq. (28) and using the decay ratethus obtained in the ratio D m /r , with the use of Eq. (24)one gets after some rearrangement D m r = (cid:18) bf (cid:19) c v w ( w /w ) f ∞ a c v w f a a . (31) Comparing this with Eq. (4) and using Eq. (8) from Eq.(27) one arrives at the expression λ | w ≫ w = Aa r DD H , (32)where A = p bf / ≈ .
111 (33)is a numerical constant. Thus, the parameter λ of thephenomenological theory in the large- w case can be cal-culated from only experimentally measurable quantities. B. Diffusion profiles of unstable I-v pairs
One consequence of the pairs instability is that theirdiffusion cannot be described by the conventional diffu-sion equation Eq. (23) because the second Fick’s law ex-presses the conservation of the diffusing particles whichis not the case with the pair diffusion. Being unstable,the I-v pairs obey instead of Eq. (23) the non-Fickiandiffusion equation suggested in Ref. 13: ∂G p ( R , t ) /∂t = D m ∇ G p ( R , t ) − rG p ( R , t ) , (34)where we introduced the pair Green’s function (GF) thatsatisfies the delta-function initial condition G p ( R , t = 0) = δ ( R ) (35)and describes the probability to find the pair at point R at time t . Integrating Eq. (34) over the space variablesit is seen that the probability to find the pair at time t anywhere in the system diminishes as e − rt , as expected.Explicit expressions for the GF of the decaying pair inthe space-time variables can be written straightforwardlyfor any dimension dG p ( R , t ) = 1(4 πD m t ) d/ exp (cid:18) − R D m t − rt (cid:19) . (36)Below we will also need the LF-transformed G p that isalso easily found from Eqs. (34) and (35) as G p ( K , z ) = 1 z + r + D m K . (37)Despite instability of the pairs, the impurity densityshould conserve irrespective of the diffusion mechanism.This is indeed the case if we take into account the im-purities from the decayed pairs that simply immobilize(become stable) and their density grows with time as G sp ( R , t ) = r Z t dt ′ G p ( R , t ′ ) . (38)Integrating this over the spacial variables is easy to checkthat the normalization of the sum G p ( R , t ) + G sp ( R , t ) (39)is equal to unity at all t .We note that the first term in Eq. (39) at large timegoes to zero so the impurity distribution is dominated bythe second term. Due to the specifics of the STM tech-nique, only the second term was observed experimentallyin Refs. 12 and 13 and only at t = ∞ in which case thefirst term vanished and the second acquired the exponen-tial asymptotic behavior in the spatial variables G sp | | R |→∞ ∝ exp (cid:0) − | R | /λ (cid:1) . (40)It can be shown that this is a universal behavior at t = ∞ in all dimensions. To see this we first notice thatthe Laplace transform of a function at z = 0 is just theintegral of the function over the time variable from zeroto infinity. Thus, using Eqs. (37) and (38) one gets withthe use of the inverse Fourier transform G sp ( R , ∞ ) = r (2 π ) d Z d d K e i K · R r + D m K = 1(2 π ) d Z d d K e i K · R λ K ) ≡ G P ( R ) , (41)where on the second line we introduced the kernel of thescreened Poisson equation G P which is known to havethe exponential asymptotic behavior Eq. (40) in all di-mensions.In contrast to the diffusion in surface layers where itis possible to observe individual I-v pairs and ignore therest of the impurities, in 3D diffusion all impurityatoms contribute so at t = ∞ the profile will span thewhole crystal and will be seen as just the homogeneousequilibrium distribution. Therefore, NGDPs can be ob-served only at finite t < ∞ in which case their asymptoticwill be Gaussian as can be easily illustrated in 1D geom-etry. Setting in Eqs. (36) and (38) d = 1 and R = X andtaking the integral over t one gets G sp ( X, t ) = r Z t G p ( X, t ′ ) dt ′ = − λ X s = ± se s | X | /λ erfc (cid:18) | X | √ D m t + s √ rt (cid:19) . (42)Because at x = −∞ erfc( x )=2, the behavior of this ex-pression at t = ∞ is exponential in | X | , as expected. Atfinite t , however, the behavior is Gaussian, as can be seenfrom the behavior of the the erfc function at large valuesof its argument erfc( x ) | x →∞ ∼ e − x / ( √ πx ) . Thus, strictly speaking the exponential tails can neverbe observed in 3D concentration profiles, only someexponentially-looking transient features, as will be shownin more general case of multiple-encounter diffusion inSec. V. Still, the profile can be very close to the expo-nential shape at low temperature and a sufficiently largedensity of the associated pairs in the initial state. In thiscase the existing pairs start to diffuse immediately while the immobile impurities need first to enter into associa-tion with vacancies. At low temperature the waiting timemay be quite long because of the small vacancy concen-tration so the impurity distribution due to the pair diffu-sion may advance to the stage where it will only slightlydiffer from the infinite-time exponential profile.
IV. IMPURITY DIFFUSION VIA MULTIPLEI-V ENCOUNTERS
In the previous section we discussed diffusion of indi-vidual I-v pairs. In particular, it was noted that Eqs.(38) and (42) that describe the impurity distributionsdue to the decayed pairs can be studied experimen-tally in the surface layers by means of the STM mi-croscopy which makes possible investigation of each I-vencounter individually disregarding the diffusion of otherimpurities.
In conventional experiments on diffusionin 3D bulk, however, such separation is not possible. Theimpurities in the profile are indistinguishable and thereis no way to differentiate them according to their evo-lution history. Therefore, theoretical description shouldtake into account all possible impurities: those belongingto the initial profile, associated impurities in the mobilestate, or the impurities that have already undergone oneor more I-v encounters. All these contributions shouldbe accounted for in a single diffusion profile with appro-priate weights. In Ref. 4 this problem was solved by firstfinding the solution for zero and one encounter and theniterating the distribution obtained as many times as nec-essary to describe the profile at a desired stage of theevolution. In the present paper we will essentially followthis route but making it more formally refined.Namely, we are going to use a technique of the many-body theory usually referred to as the Dyson equation(see, e. g., Ref. 51). In this approach the problem ofrepeated interactions of a particle is separated into anirreducible and a reducible parts and the repeated itera-tions of the irreducible part are simply summed up as ageometric series as G + G Σ G + G Σ G Σ G + · · · = 1 G − − Σ . (43)Here the products are either convolutions in the space-time variables or the usual algebraic products of the LF-transformed quantities; G is the GF of a free particleand Σ is the irreducible interaction part that cannotbe represented as two interactions separated by the freepropagation, as, e. g., in the third term on the left handside (l.h.s.). The practical observation that makes thisapproach useful is that the combinatorial problem of offinding the sum of all contributions that include the freepropagation and a single or multiple I-v encounters asrepresented in compact form on the r.h.s. of Eq. (43)can be fully recovered from only the first two terms onthe l.h.s. with the first term being known exactly. Thisobservation was successfully applied to the problem ofvacancy-mediated diffusion in Refs. 52 and 53. The au-thors effectively derived expressions for Σ in Eq. (43)for the cases of the self-diffusion and in a two-frequencymodel and we are going to apply this approach to the5FM. We will call the irreducible part Σ the diffusionkernel and express it through another quantity D thatwill be called the diffusivity as followsΣ( K , z ) = −D ( K , z ) K . (44)The free GF G is easily found from the observation thatwithout interaction with the vacancies the impurity isimmobile and remains in its initial position: G ( R , t ) = G ( R , t = 0) = δ ( R ). With the use of the LF transformone easily finds G ( z ) = Z ∞ dt e − zt Z d d R e − i K · R δ ( R ) = 1 z . (45)The impurity GF with all I-v encounters being taken intoaccount is obtained by substituting the last two equationsEqs. (44) and (45) into Eq. (43): G ( K , z ) = 1 z + D ( K , z ) K ≈ z − z D ( K , z ) K . (46)The diffusion constant is found from the diffusivity as D = D ( K = , z = 0) . (47)But our main interest is in a nontrivial dependence of D on its arguments K and z because if the diffusivitywere independent of K and z the GF in Eq. (46) whentransformed to the time and space variables would bestrictly Gaussian without any traces of the NGDPs weare interested in.Thus, our goal is to find explicitly the last term in Eq.(46) and with its use to recover the complete impurityGF. To to fulfill this goal we will use Eqs. (3) and (4)from Ref. 4 that in the GF notation read ∂G m /∂t = D m ∇ G m − rG m + gG s ∂ ( G s + G m ) /∂t = ∂G/∂t = D m ∇ G m , (48)where the total impurity GF G = G m + G s is separatedinto the mobile ( G m ) and immobile ( G s ) parts with thesubscript “ s ” standing for “static” because in contrast toRefs. 4 and 5 in the 5FM the impurity is always in thesubstitutional position. Parameter g in Eq. (48) is therate of transition of the impurity from the static to themobile state which in the 5FM is the rate of I-v associ-ation. Because in the association participates a vacancy,the rate g should be proportional to the vacancy concen-tration and thus is of order O ( c v ). Its calculation will bediscussed below.The LF-transformed set of Eqs. (48) is( z + r + D m K ) G m ( K , z ) − gG s ( K , z ) = 12 c NN ( z + D m K ) G m ( K , z ) + zG s ( K , z ) = 1 , (49) where we assumed the initial conditions G m ( R , t = 0) = 12 c NN δ ( R ) G s ( R , t = 0) = (1 − c NN ) δ ( R ) G ( R , t = 0) = δ ( R ) . (50)Here c NN is the vacancy concentration at the NN sitesof the impurity which can be different from c v due tothe I-v interaction or because of the way the initial statewas prepared, though we will assume that it is of thesame order of magnitude as c v ; 12 c NN in Eqs. (50) isthe density of associated I-v pairs in the initial state inthe FCC lattice which coordination number is 12 whichimplicitly assumes that all sites around the vacancy areassumed to be equivalent.It is important to note that in applying the Fouriertransform we assume that the crystal is translationallyinvariant which in particular means that the rate of as-sociation g is a position-independent constant. But be-cause g depends on the vacancy concentration, this maynot be the case in experiments where non-equilibrium va-cancy concentration may acquire inhomogeneity becauseof their influx from the surface or deposited layers. Such cases cannot be treated by Eqs. (49). It is to beunderstood that we are considering the dilute systemsand c v should be constant far from the impurities in thehost bulk. In the vicinity of impurity it can be differ-ent from its bulk value due to the I-v interaction bothat and out of thermal equilibrium. In the latter case inthe Smoluchowski picture of the vacancy capture on theimpurity NN sites a non-constant vacancy diffusion pro-file forms near the impurity because the NN sites serveas the sinks. Because, as we pointed out in Sec. II, the 5FM ade-quately describes the I-v interaction only to first order inthe vacancy concentration, we will solve the system Eqs.(49) only to this order by first finding the mobile GF G m = 12 c NN z + gz ( z + r + D m K ) . (51)With known G m the total GF G = G s + G m can be founddirectly from the second of Eqs. (49) as G = 1 z − z (12 c NN z + g ) D m K z + r + D m K . (52)As is seen, the approach Ref. 4 turned out to be similarto that of Refs. 52 and 53 by giving only the zeroth andthe first order terms of the expansion in Eq. (43) so toobtain the full impurity GF the Dyson equation will haveto be used. Comparing Eq. (52) with Eq. (46) we arriveat the expression for the K and z -dependent diffusivity D ( K , z ) = 12 c NN z + gz + r + D m K D m . (53)The two contribution to the second term in Eq. (52) de-scribe somewhat different diffusion scenarios so let us dis-cuss them separately. The term proportional to c NN de-scribes the associated I-v pairs that are present in theinitial profile. Their density 12 c NN can be arbitrary de-pending on the way the profile was prepared. From Eqs.(53) and (47) one can see that this term does not con-tribute to the diffusion constant which is natural becausean arbitrary initial condition cannot influence the quan-tity corresponding to thermal equilibrium. The physicalmeaning of this term becomes transparent after its rear-rangement into three contributions G | c NN = − c NN D m K z ( z + r + D m K ) = − c NN z + 12 c NN z + r + D m K + 12 c NN rz ( z + r + D m K ) . (54)From Eq. (37) it can be seen that apart from the factor12 c NN , the first term on the second line in Eq. (54) isthe LF transformed G p Eq. (37) while the second term isEq. (37) multiplied by r/z . But multiplication by 1 /z ofthe Laplace transform of a function corresponds to theLaplace transform of the integral of this function over t . Thus, the second term on the second line of Eq. (54)corresponds to LF-transformed G sp from Eq. (38), so thesum of the two terms is Eq. (39) multiplied by 12 c NN .Thus, these terms describe the diffusion profile of thepairs of density 12 c NN that were present in the initialstate. The negative term on the first line in Eq. (54)simply accounts for the fact that the associated impuri-ties were taken from the initial delta-function profile (seeEqs. (50) and (45)).In contrast to the term that describes diffusion of theI-v pairs which already exist in the initial state, the termproportional to g in Eq. (52) describes the diffusion of ini-tially immobile impurities. To enter into the mobile statethe impurities need first to be associated with a vacancy.This process is limited by the low vacancy concentrationand so the ensuing diffusion is much slower than the pairdiffusion. The time scales of the two diffusion modes aredefined by the lifetime of the pairs τ decay = 1 /r and bythe characteristic time of the I-v association τ assn. = 1 /g .To compare their relative values we first note that sub-stitution of diffusivity Eq. (53) into Eq. (47) gives D = gr D m (55)which leads to the relation τ assn. τ decay = D m D ≫ , (56)where the last inequality follows from Eq. (25) and thediscussion that follows it. At low temperatures this ratiocan be large. In the LaAl system at 50 ◦ C we estimatedit in Eq. (26) as amounting to two orders of magnitude.The main reason for this is that while the decay rate r is of zeroth order in the vacancy concentration ( O ( c v ) = O (1)), the association rate g is of O ( c v ), as can be seenfrom the equation g = rD/D m = 84 c v w p ∞ (57) obtained from Eqs. (55), (25), and (28).The expression for the rate of I-v association Eq. (57)was obtained in the framework of the 5FM while in Ref.4 a Smoluchowski-type formula was used (see their Eq.(12)). To compare the two approaches let us first considerEq. (57) in the case w = w because the Smoluchowskiformula Φ = 4 πR c ρ D v (58)describes the flux of vacancies of density ρ that are be-ing caught by the sphere of radius R c and depends onthe vacancy diffusion constant D v . The latter is definedin a continuous homogeneous medium where all diffu-sion steps are equivalent which in the lattice case meansthat w cannot be different from w . Thus, substituting p ∞ ( w = w ) = 1 / (1 + b ) from Eq. (29) in Eq. (57) weget g ( w = w ) ≃ . c v w . (59)To find the value from the Smoluchowski formula we notethat the per volume vacancy density is ρ = 4 c v /a be-cause there is four sites in the cubic cell in the FCClattice, the first CS radius R c = a/ √
2, and D v = w a .Substituted into Eq. (58) this gives g S ≈ . c v w (60)in excellent agreement with the 5FM value Eq. (59).As our analysis of the database for aluminum host hasrevealed, the case of large w , is particularly interestingfor experimental purposes. To apply the Smoluchowskiequation to this case we remind that as w → ∞ allvacancies that arrive at the sites belonging to classes 4–13 in Fig. 2 immediately form bound I-v pair (see Sec.III). This means that the capture radius R c is effectivelyshifted toward a larger value that can be assessed byaveraging the distances from the impurity to all sites inthese classes. Elementary calculation gives R c ≈ . a and from Eq. (58) g S | w ≫ w ≈ c v w . (61)which is also very close to the 5FM value ∼ c v w thatcan be obtained from Eq. (57) with the use of Eq. (30).Thus, our Eq. (57) agrees with the formula suggestedin Ref. 4 in two cases where the Smoluchowski formulais applicable but in addition covers the cases of arbitraryvalues of w . But more important to us is that Eq. (57)for large w can be cast in the form g | w ≫ w = (84 /bf ) c v w f = D H / ( Aa ) , (62)where A is given by Eq. (33). Thus, Eqs. (32) and (62)allow us to express two parameters of the phenomenolog-ical theory only in terms of experimentally measurablequantities.0 V. DIFFUSION PROFILES
In the GF approach the diffusion profiles are obtainedby convolution of the initial profile with the impurity GF.The latter is obtained in our approach by first substitut-ing the diffusivity from Eq. (53) into Eq. (46) and thentaking the inverse LF transform to find the GF in thespace and time variables. Before proceeding with con-crete implementation of this procedure we have to agreeon the value of the density of the I-v pairs in the initialprofile that was estimated to be equal to 12 c NN . Thisestimate presumes that the vacancy can be found at dif-ferent NN sites of the impurity with equal probability.That may not be the case if the vacancies are introducedin the initial profile by means of a non-equilibrium tech-nique that causes non-isotropic distribution of the vacan-cies on the NN positions of the impurities. Such casesare beyond the scope of our approach which is restricted,as we pointed out in Sec. IV, to stationary and homo-geneous distributions of I-v pairs, though not necessarilycorresponding to thermal equilibrium. The latter, how-ever, is a natural choice, so all estimates will be done forthis case. In particular, the equilibrium density of theI-v pairs can be found from the expression12 c ( eq ) NN ≃ c v e E b /k B T . (63)With E b ≈ .
47 (see Table I) the associated impuritieswill constitute at 50 ◦ C about 1% of their total num-ber which is a small but detectable quantity. However,to observe the NGDP behavior which takes place atshort distances one needs to keep initial profiles maxi-mally sharp.
But for establishment of the equilibriummany association-dissociation events must occur accom-panied by impurity diffusion with ensuing profile smear-ing. Therefore, we will assume that the initial distri-bution was prepared at a temperature so low that cor-responding c ( eq ) NN is negligible and that the preparationtechnique does not introduce excess vacancies. So in ourcalculations below we for simplicity will neglect the termsproportional to c NN in Eq. (52). In case of necessity theycan be taken into account along the lines of derivationpresented below. Another reason for omission of theseterms is that this reduces the problem to the case stud-ied in Ref. 4, thus facilitating comparison between thetwo approaches.Under approximation c NN ≃ G ( K , z ) = (cid:18) z + gD m K z + r + D m K (cid:19) − . (64)This expression can be cast into the form convenient forthe inverse Laplace transform and for assessment of therelative magnitude of different contributions: G ( K , z ) = 1 z − z + z z z − z (cid:18) z − z − z − z (cid:19) (65) where z , = − r + D m K ± r ( r + D m K ) − gD m K (66)or to the leading order in g = O ( c v ) z ≃ − gD m K r + D m K = Σ( K ,
0) (67) z ≃ − ( r + D m K ) . (68)where the diffusion kernel isΣ( K ,
0) = − gD m K r + D m K = g (cid:18)
11 + ( λ K ) − (cid:19) . (69)Neither the kernel nor the diffusivity now do not de-pend on z and in equations below this argument will bedropped.The Laplace transform in Eq. (65) reduces to the cal-culation of pole residues: G ( K , t ) = e tz + z z z − z (cid:0) e tz − e tz (cid:1) ≃ e t Σ( K ) + D ( K ) K (cid:16) e tz − e t Σ( K ) (cid:17) . (70)As can be seen, in real space the second term on the sec-ond line would integrate to zero while the first one tounity because the spacial integration corresponds to theFourier component K = . This separates Eq. (70) intocontributions of different order in c v as follows. Follow-ing Cowern et al., let us consider the long-time dif-fusion when the number of impurities from the initialprofile that experienced one or more encounters with thevacancies is comparable to the number of all impuritiesin the profile, i. e., is of order O (1). This means that gt = O (1) which in its turn means t = O (1 /c v ). Nowbecause D = O ( c v ), the products t D in the first and inthe last exponential functions are of order unity. Butthe factor D before the second term makes the contribu-tion due to the last exponential function negligible if weare interested only in O (1) terms. The only potentiallyproblematic term is the first exponential function in theparentheses that contains t = O ( c − v ) without compen-sating O ( c v ) factor because z = O (1) as can be seenfrom Eq. (68). But as is easy to see, e tz ≤ e − rt ≪ r = O (1) while t = O ( c − v ). Thus, the secondterm in Eq. (70) is much smaller than the first term andcan be neglected if we agree to neglect O ( c v ) contribu-tions.Thus, the leading O (1) terms in Eq. (70) can bereduced to the following Fourier-transformed diffusionequation ∂C ( K , t ) ∂t = Σ( K ) C ( K , t ) , (72)where C is the Fourier transform of impurity concen-tration. In real space the diffusion equation Eq. (72)1turns out to be not a differential equation but an integro-differential one ∂C ( R , t ) ∂t = Z Σ( R − R ) C ( R , t ) d R (73)with the diffusion kernelΣ( R − R ) = gG P ( R − R ) − gδ ( R − R ) . (74)Here G P is the GF of the screened Poisson’s equationfrom Eq. (41) and thus describes the limiting ( t → ∞ )profile of the pair diffusion. This makes transparent thephysical meaning of Eq. (73). The impurity at point R is picked up by a vacancy with the rate g (the second termon the r.h.s.) and via the pair diffusion is redistributedwith the probability density G P (the first term on ther.h.s.). The time dependence of this process is ignoredbecause Eq. (73) describes diffusion on the time scale O ( c − v ) which is much larger than the O ( r − ) = O (1)scale of the pair diffusion. This is exactly the physicsstudied in Refs. 4 and 5. From Eq. (73) is easy to un-derstand the non-Fickian character of the pair-mediateddiffusion discussed in Refs. 4 and 11. Imagine a hostwith inhomogeneous distribution of impurities in it and abounded region within that is completely devoid of them.Despite this, the rate of growth of the impurity concen-tration inside the region will be everywhere positive ac-cording to Eq. (73) because of the spatially extended dif-fusion kernel Eq. (74) that is able to displace impuritiesat finite distances. This starkly contrasts with the localcurrent picture underlying the Fickian diffusion.3D diffusion equation Eq. (73) makes possible numeri-cal simulation of the diffusion profiles in any geometry. Inlarge systems the solution should presumably be soughtvia its direct numerical integration. In 3D the Poissonkernel is singular (it coincides with the screened Coulombor the Yukawa potential) but efficient techniques of deal-ing with it were proposed in Ref. 49. In systems of moder-ate sizes a convenient method provides the Fourier trans-form. In 1D case the solution reads C ( X, t ) = 12 π Z ∞−∞ e t Σ( K ) C ( K ) dK, (75)where C is the Fourier transform of the initial profile.In this way were calculated NGDPs presented in Fig. 4where C where chosen to be Gaussian to easier visual-ize the NGDPs caused by the pair-mediated diffusion (inconventional diffusion an initially Gaussian-shaped pro-file remain Gaussian at all times). As shown in AppendixA, the 1D profiles calculated within our approach coin-cide with those of Refs. 4 and 5 so in Fig. 4 one can seethe exponential tails in the diffusion profiles similar tothose found by Cowern et al. The tails, however, cannotextend on arbitrarily long distances because the diffusionkernel in Eq. (69) at small Fourier momenta behaves as ∼ K so from the inverse Fourier transform of the kind ofEq. (75) but for arbitrary dimension it is easy to see thatat finite times the large- | R | asymptotic will be Gaussian,similar to the case of diffusion of individual pair discussedin Sec. III B. C ( X , t )( a r b . un i t s ) X ( µ m) FIG. 4. 1D profiles for the diffusion of lanthanum impurityin aluminum at temperature 50 ◦ C starting from the initialGaussian distribution of with 0.5 µm (dashed line). Thicksolid line: the profile calculated according to Eq. (75) for t = 5 hours; thin solid line: t = 50 hours; dashed-dottedlines: the profiles at 5 and 50 hours as predicted by the con-ventional diffusion equation; the time intervals were chosento correspond to gt ≃ ≃
10, respectively.
VI. CONCLUSION
In the present paper a theory of the vacancy-mediateddiffusion in the case of strong I-v binding has been de-veloped. It has been shown that tightly bound I-v pairsprovide the mobile state of impurities that underly theNGDPs similar to those observed in dopant diffusion insemiconductors and in the copper surface layers.
By unifying the phenomenological theory of Cowern etal. with the 5FM of the vacancy-mediated diffusion inFCC hosts it has been possible to calculate numericalvalues of the parameters of the phenomenological theoryon the basis of the available data on the parameters ofthe 5FM. This has made possible identificationof the impurity-host systems suitable for the observationof the NGDPs as well as their explicit simulation in LaAlsystem where the phenomenon is expected to be the mostpronounced among the solutes in aluminum host. Be-cause the NGDPs are universal, all impurities listed inTable II should exhibit the same profiles as shown in Fig.4 but at shorter length scales. The latter can be enlargedby lowering the temperature but the time of the obser-vation will have to be extended correspondingly.Apart from the calculation of the parameters of thephenomenological model in the framework of the 5FM,the approach of Cowern et al. has been extended intwo respects. First, it has been shown that the diffusion2mode studied in Refs. 12 and 13 in 2D can contribute to3D NGDPs in cases when the initial state already con-tains associated I-v pairs. Because their diffusion startsimmediately, they introduce impurity diffusion on muchshorter time scale than that of Ref. 4. For example, inLaAl system at room temperature (20 ◦ C) the parame-ter λ ≃ µm . The conventional diffusion will cower thisdistance in over a month while the pair diffusion in aboutthree hours. Accounting for this mode of diffusion maybe of practical importance in assessment of the longevityof microelectronic devices. Though the concentration ofpreexisting pairs is usually expected to be small, in someof the projected devices the functional elements will con-sist of only one atom so the estimates of the longevitythat neglect the fast diffusion of the small number of con-taminating atoms that are associated with the vacanciesintroduced during the deposition process may lead to se-rious errors.The second extension of the theory of Ref. 4 has beenachieved through its blending with the Dyson equation.This resulted in a non-Fickian integro-differential diffu-sion equation describing the pair-propagated impuritydiffusion in arbitrary geometry that can be used in sim-ulations of NGDPs in any elemental hosts.Special attention in the paper has been devoted to thesystems with large w /w ratio for two reasons. First,the ratio turned out to be large for all impurities in alu-minum host with the largest values of λ , i. e., in thesystems that should be the most appropriate for exper-imental study of NGDPs. Secondly, and more impor-tantly, this case makes possible accurate quantitative pre-dictions about the phenomenon. A serious problem ofthe microscopic diffusion theory is that both experimen-tal definitions and first-principle calculations of variousactivation energies contain errors of order O (0 . The jump frequencies and diffusionconstants depend on the energies via the Arrhenius lawand at temperatures in a few hundred Kelvins may beorders of magnitude off from their true values. But ithas been shown in the present paper that when w /w is large the phenomenological parameters g and λ can becalculated from only two quantities: the impurity diffu-sion constant and that of the host self-diffusion. Both canbe measured at the experimental temperature indepen-dently and the two parameters calculated on their basiscan subsequently be used in the profile simulations. Thesimulated NGDPs should agree with experimental onesquantitatively, provided the 5FM is an adequate modelfor the system under study. Significant discrepancies willmean that the canonical 5FM is too simplistic for the caseunder consideration and needs to be improved along thelines suggested by first-principles calculations and phys-ical considerations. This conclusion relies on theassumption that the errors in in the frequencies w and w are not too large to reduce our estimates of w /w ra-tios for the systems listed in Table I more than 4–6 ordersof magnitude. Unfortunately, at present this possibilitycannot be completely excluded. The important question that has not be adequatelyaddressed in the present paper concerns the reliabilityand the accuracy of the developed theory which has beensubstantiated mainly by qualitative arguments and phe-nomenological approaches. This difficult question will beaddressed in a separate paper where a rigorous treatmentof the 5FM in the general case of the I-v interaction ofarbitrary strength will be presented. It will be shownthat in the limit of strong I-v attraction the results ofthe present paper are in excellent agreement with therigorous solution. ACKNOWLEDGMENTS
I would like to express my gratitude to Hugues Dreyss´efor encouragement.
Appendix A: Comparison with NGDP of Ref. 4
In slightly modified notation, the 1D diffusion profilewith initial delta-function distribution was shown to bedescribed by the series given by Eqs. (5)-(10) of Ref. 4 as C ( ξ, θ ) = λ − ∞ X n =0 P n ( θ ) φ n ( ξ, , (A1)where θ = gt, ξ = X/λ, (A2) P n ( θ ) = ( θ n /n !) e − θ , (A3) φ n =0 ( ξ,
1) = δ ( ξ ) , (A4)and φ n> ( ξ, µ ) = e −√ µ | ξ | (2 √ µ ) n − n − X k =0 k k ! (cid:18) n − − kn − (cid:19) ( | ξ |√ µ ) k . (A5)Here we introduced the parameter µ to facilitate theproof that the profile Eq. (A1) from Ref. 4 coincideswith the 1D profile from Eq. (75) with C ( K ) = 1 forthe delta-function initial profile: C ( X, t ) = 12 π Z ∞−∞ dKe iXK e − gt exp (cid:18) gt λK ) (cid:19) = 1 λ ∞ X n =0 P n ( θ ) 12 π Z ∞−∞ e iξu du (1 + u ) n . (A6)Here the last exponential on the first line has been ex-panded in the Tailor series so by comparison with Eq.(A1) we conclude that the inverse Fourier transforms onthe second line should be equal to φ n> ( ξ, φ n ( ξ, µ ) = 12 π Z ∞−∞ e iξu du ( µ + u ) n (A7)3and note that if φ n =1 ( ξ, µ ) is known, other integrals canbe computed recursively as φ n +1 ( ξ, µ ) = − n ddµ φ n ( ξ, µ ) . (A8)Thus, we only need to show that φ n ( ξ, µ ) in Eq. (A5)satisfy the recursion. To this end we first note that withthe exponential factor being common to all terms in allfunctions, the equality in Eq. (A8) will hold if it will bevalid for every power of | ξ | k under the summation sign. Let us consider one such term in Eq. (A5) φ ( k ) n = e −√ µ | ξ | n − ( n − k k ! (2 n − − k )!( n − − k )! | ξ | k µ ( k +1) / − n . (A9)When substituted in Eq. (A8) it will contribute to | ξ | k term in φ ( k ) n +1 through the derivative of its last factorwith respect to µ . The only other contribution from φ n contributing into | ξ | k term in φ n +1 is φ ( k − n differenti-ated with respect to µ in the exponential function. Itis straightforward to check that these two contributionslead to the term φ ( k ) n +1 as in Eq. (A9) only with n + 1instead of n , as required. P. A. Packan, MRS Bulletin , 1821 (2000). J. C. Ho, R. Yerushalmi, Z. A. Jacobson, Z. Fan, R. L.Alley, and A. Javey, Nature Materials , 62 (2008). P. M. Fahey, P. B. Griffin, and J. D. Plummer, Rev. Mod.Phys. , 289 (1989). N. E. B. Cowern, K. T. F. Janssen, G. F. A. van de Walle,and D. J. Gravesteijn, Phys. Rev. Lett. , 2434 (1990). N. E. B. Cowern, G. F. A. van de Walle, D. J. Gravesteijn,and C. J. Vriezema, Phys. Rev. Lett. , 212 (1991). C. S. Nichols, C. G. Van de Walle, and S. T. Pantelides,Phys. Rev. Lett. , 1049 (1989). M. Y. L. Jung, R. Gunawan, R. D. Braatz, and E. G.Seebauer, AIChE Journal , 3248 (2004). K. Chen, R. Vaidyanathan, E. G. Seebauer, and R. D.Braatz, J. Appl. Phys. , 026101 (2010). R. Vaidyanathan, M. Y. L. Jung, R. D. Braatz, and E. G.Seebauer, AIChE Journal , 366 (2006). R. Kube, H. Bracht, E. H¨uger, H. Schmidt, J. L. Hansen,A. N. Larsen, J. W. Ager, E. E. Haller, T. Geue, andJ. Stahn, Phys. Rev. B , 085206 (2013). S. Mirabella, D. De Salvador, E. Napolitani, E. Bruno,and F. Priolo, J. Appl. Phys. , 031101 (2013). R. van Gastel, E. Somfai, W. van Saarloos, and J. W. M.Frenken, Nature , 665 (2000). R. van Gastel, E. Somfai, S. B. van Albada, W. van Saar-loos, and J. W. M. Frenken, Phys. Rev. Lett. , 1562(2001). R. van Gastel, R. Van Moere, H. J. W. Zandvliet, andB. Poelsema, Surface Science , 1956 (2011). M. L. Grant, B. S. Swartzentruber, N. C. Bartelt, andJ. B. Hannon, Phys. Rev. Lett. , 4588 (2001). M. J. A. M. Brummelhuis and H. J. Hilhorst, J. Stat. Phys. , 249 (1988). Z. Toroczkai, Int. J. Mod. Phys. B , 3343 (1997). A. Chroneos, H. Bracht, R. W. Grimes, and B. P. Uberu-aga, Applied Physics Letters , 172103 (2008). A. B. Lidiard, Phil. Mag. , 1218 (1955). A. L. Claire, J. Nucl. Mater. , 70 (1978). J. Manning,
Diffusion Kinetics for Atoms in Crystals (VanNostrand, Princeton, N. J., 1968). J. Philibert,
Atom Movements (Les ´Editions de Physique,Les Ulis, 1991). S. Mantl, W. Petry, K. Schroeder, and G. Vogl, Phys. Rev.B , 5313 (1983). U. Klemradt, B. Drittler, T. Hoshino, R. Zeller, P. H. Ded-erichs, and N. Stefanou, Phys. Rev. B , 9487 (1991). M. Mantina, S. L. Shang, Y. Wang, L. Q. Chen, and Z. K.Liu, Phys. Rev. B , 184111 (2009). C. Wolverton, Acta Materialia , 5867 (2007). D. Simonovic and M. H. F. Sluiter, Phys. Rev. B ,054304 (2009). M. Mantina, Y. Wang, L. Chen, Z. Liu, and C. Wolverton,Acta Materialia , 4102 (2009). C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer,G. Kresse, A. Janotti, and C. G. Van de Walle, Rev. Mod.Phys. , 253 (2014). H. Wu, T. Mayeshiba, and D. Morgan, Scientific Data ,160054 (2016). H. Wu, T. Mayeshiba, and D. Morgan, (2017),10.6084/m9.figshare.1546772.v8. T. Angsten, T. Mayeshiba, H. Wu, and D. Morgan, NewJ. Phys. , 015018 (2014). M. Krivoglaz, JETP , 1273 (1961). M. A. Krivoglaz and S. P. Repetskiy, Fiz. Met. Met-alloved. , 899 (1971), [Phys. Met. Metallogr. (USSR)32, 1 (1971)]. O. B´enichou and G. Oshanin, Phys. Rev. E , 020103(2001). J. R. Manning, Phys. Rev. , 2169 (1962). C. Tuijn, H. Bakker, and G. Neumann, J. Phys. Condens.Matter , 4801 (1992). J. L. Bocquet, Phil. Mag. , 3603 (2014). C. Zacherl, S.-L. Shang, and Z.-K. Liu, (2012). K. Carling, G. Wahnstr¨om, T. R. Mattsson, A. E. Matts-son, N. Sandberg, and G. Grimvall, Phys. Rev. Lett. ,3862 (2000). J. R. Manning, Phys. Rev. , A 1758 (1964). M. Koiwa and S. Ishioka, J. Stat. Phys. , 477 (1983). Y. Du, Y. Chang, B. Huang, W. Gong, Z. Jin, H. Xu,Z. Yuan, Y. Liu, Y. He, and F.-Y. Xie, Mater. Sci. Engin.:A , 140 (2003). A. B. Bortz, M. H. Kalos, and J. L. Lebowitz, J. Comput.Phys. , 10 (1975). N. G. van Kampen,
Stochastic Processes in Physics andChemistry (North-Holland, NY, 1981). R. M. Ziff, S. N. Majumdar, and A. Comtet, J. Chem.Phys. , 204104 (2009). V. I. Tokar, (2005), arXiv:cond-mat/0505019. V. I. Tokar, (in preparation). A. Cerioni, L. Genovese, A. Mirone, and V. A. Sole, J.Chem. Phys. , 134108 (2012). E. W. Weisstein, “Erfc.” http://mathworld.wolfram.com/Erfc.html (visited on 12/13/2017). D. Thouless,
The Quantum Mechanics of Many-Body Sys-tems: Second Edition , Dover Books on Physics (Dover Publications, 2014). R. A. Tahir-Kheli and R. J. Ellott, J. Phys. C , L445(1982). R. A. Tahir-Kheli and R. J. Elliott, Phys. Rev. B , 844(1983). C. S. Nichols, C. G. Van de Walle, and S. T. Pantelides,Phys. Rev. Lett. , 1049 (1989). E. W. Weisstein, “Laplace Transform.” http://mathworld.wolfram.com/LaplaceTransform.html . B. E. Kane, Nature393