Impurity Green’s function in the five frequency model of diffusion in the FCC host: General case and the limit of strong impurity-vacancy binding

Impurity Green’s function in the five frequency model of diffusion in the FCC host: General case and the limit of strong impurity-vacancy binding

V. I. Tokar Université de Strasbourg, CNRS, IPCMS, UMR 7504, F-67000 Strasbourg, France
July 13, 2019

The impurity Green’s function exact to the first order in the vacancy concentration has been calculated in the framework of the five frequency model (5FM). The solution in terms of determinant ratios has been obtained with the use of the Cramer’s rule. The determinant sizes varied from 54 for the most general case to three for the four frequency model. Both analytical and numerical techniques has been used to analyze the solution. Special attention has been devoted to the case of strong impurity-vacancy (I-v) binding in order to substantiate the picture of diffusion via bound I-v pairs developed earlier in a phenomenological approach. Complete agreement with the phenomenological theory has been established thus providing its rigorous justification. The solution has been also applied to the calculation of the diffusional broadening of the Mössbauer resonance in FeAl system and good agreement with available experimental data and calculations within the encounter model has been found. The decay of the density of the nearest neighbor I-v pairs has been discussed in detail and suggested to be used as an additional constraint on the parameters of the 5FM. The possibility of experimental observation of the decay with the use of the positron annihilation technique has been pointed out.

I Introduction

In Ref. Tokar, 2018 (in the following referred to as I) the impurity Green’s function (GF) of the five-frequency model (5FM) of the vacancy-mediated diffusionLidiard (1955); Claire (1978) was calculated in the limit of a strong impurity-vacancy (I-v) attraction. The solution was based on the phenomenological approaches of Refs. Cowern et al., 1990; van Gastel et al., 2001 where the bound I-v pairs were treated as quasiparticles which diffusion was shown to be responsible for the non-Gaussian diffusion profiles (NGDPs) and the non-Fickian diffusion experimentally observed in the bulk diffusion in siliconCowern et al. (1990) and in the Cu(001) surface layer.van Gastel et al. (2000) Because of the instability of the bound pairs, the diffusion proceeds via repeated I-v associations and subsequent decays with the encounters separated by long periods of immobility because of the low vacancy concentration that usually does not exceed 0.1 at.%̃ even at the melting point.Angsten et al. (2014); Carling et al. (2000) After large number of encounters the impurity density distribution acquires the conventional Gaussian shape but at the early stage of diffusion it exhibits the exponential-looking tails characteristic of a defect-mediated diffusion.Brummelhuis and Hilhorst (1988); Cowern et al. (1990); Toroczkai (1997); van Gastel et al. (2000)

In I the phenomenological approach was implemented within the 5FM with the model parameters taken from the first-principles calculations of Refs. Wu et al., 2016, 2017 and the universal NGDPsCowern et al. (1990, 1991) were simulated. However, in Refs. Brummelhuis and Hilhorst, 1988; Toroczkai, 1997; Grant et al., 2001; van Gastel et al., 2011 it was shown both theoretically and experimentally that for the existence of exponential profiles the I-v attraction is not indispensable and may occur even under the I-v repulsion.van Gastel et al. (2011); Grant et al. (2001) Furthermore, the picture of bound I-v pairs diffusing at distances much larger than the lattice constant does not accord well with the numerical simulations of Ref. Mantl et al., 1983 where despite rather strong I-v binding (0.29 eV) the I-v encounter amounted on average to only two exchanges between the impurity and the vacancy before their separation. This is only 50% more than the tracer-vacancy exchanges in the self-diffusion with no I-v attraction. So the question arises about the uniqueness of interpretation of NGDPs in terms of the pair diffusion and of interrelation of this mechanism with the alternatives that do not require I-v attraction, as well as of the ways of making distinction between them.

The aim of the present paper is to justify the pair diffusion picture suggested in I within a rigorous master equation and Green’s function formalism.Koiwa and Ishioka (1983); Tokar (2005) Our approach will be based on generalization of the technique developed in Refs. Tahir-Kheli and Ellott, 1982; Tahir-Kheli and Elliott, 1983 for the cases of self-diffusion and of the diffusion in a two frequency model. The method of solution that will be developed in the present paper can also be applied to diffusion in any elemental crystalline hosts and to more complex models with further than nearest neighbor (NN) vacancy jumps and with arbitrary but finite I-v interaction range. Such extensions has been suggested in literature on the basis of physical arguments and the first-principles calculations.Claire (1978); Wu et al. (2016); Bocquet (2014); Tuijn et al. (1992)

The text of the paper is organized as follows. In the next section an infinite system of rate equations for the impurity GF accurate to the first order in the vacancy concentration will be derived; in Sec. III it will be reduced to a finite system of 54 equations amenable to solution by Cramer’s rule. In Sec. IV the solution will be used for illustrative calculation of the diffusional broadening of the Mössbauer resonanceMantl et al. (1983); Steinmetz et al. (1986) and for discussion in this context of the encounter model (EM).Wolf (1977); Bender and Schroeder (1979); Vogl (1996); Mantl et al. (1983) In Sec. V the finite system obtained in Sec. III will be farther reduces to the system of 13 equations which are sufficient for the study of macroscopic diffusion. In Sec. VI the diffusion constant and the parameters of the tightly bound I-v pairs will be studied numerically to confirm the analytic expressions derived in the phenomenological approach. In Sec. VII the complications that arise in the case of I-v interaction of arbitrary strength will be discussed and a simpler four frequency model (4FM) that makes possible explicit analytic solution for impurity GF will be introduced. Dissociation of NN pairs which may be present in the initial state will be discussed in detail and a possibility of experimental measurement of the evolution will be suggested. It will be shown that the data on the NN density evolution can provide additional information about the values of the 5FM parameters. In final Sec. VIII the main results will be briefly summarized.

Ii Rate equations for the impurity GF

In the 5FM all interactions in the system comprising the host crystal, the solute impurities, and the vacancies are reduced to interactions within one I-v pair.Lidiard (1955); Claire (1978); Philibert (1991); Tokar (2018) Hence, the physics that the model can describe is that of diffusion in a dilute alloy where each impurity can be treated independently, so theoretical consideration can be restricted to only one of them. Because during one I-v encounter the impurity is displaced only on a microscopic distance, for macroscopic diffusion a large number of repeated encounters is necessary.Wolf (1977) Therefore, the “minimal model” that will be studied below is a host crystal consisting of sites containing one substitutional impurity and vacancies distributed within the crystal with concentration . The large number of vacancies makes the problem a truly many-body one but the small vacancy concentration makes the solution of the first order in in most cases to be sufficient. Besides, the interactions of impurity with two and more vacancies would be very different from the simple superposition of independent I-v interactions and would require a large number of additional parameters to be introduced, such as the values of multi-vacancy interactions. This would strongly diminish the usefulness of the modelMantl et al. (1983) so in the absence of adequate description of the impurity interacting with more than one vacancy the physically meaningful solution should not go beyond the pair I-v interactions. This reduces the task to essentially a two-body problem that can be solved exactly. Such a solution will be presented below and from a physical point of view it may be considered as the best possible one because it cannot be improved further without redefinition of the model.

Because 5FM is a stochastic model, it can be exactly described by a suitable master equation (ME) (see Ref. van Kampen, 1981 and Appendix A). Its direct solution, however, usually is efficient only in the case of a few-body system, like the bound two-particle I-v pair treated in I with the use of ME. In the case of the infinite number of vacancies a straightforward use of ME would require dealing with the coordinates of all vacancies in the system which is superfluous and practically unfeasible. An appropriate way of proceeding is to average over positions of all vacancies that do not participate in the I-v encounter. This can be achieved within the rate equations (REs) approach which is, in fact, a method of solution of the ME, as explained in Appendix A.

To derive the necessary REs let us define the impurity GF as the probability density of finding the impurity at time at the lattice site described by vector with the integer coordinates if at time it was located at site


This choice of the initial position was made for convenience because GF for arbitrary initial position can be trivially obtained by the substitution .

To express GF as a statistical average let us introduce the occupation number which is equal to unity when site is occupied by the impurity and is zero otherwise. In this notation


where the angular brackets denote the average over the time-dependent nonequilibrium statistical ensemble defined in Eq. (97). Thus, besides vacancies the ensemble contains the inhomogeneous impurity density distribution . The evolution continues until the impurity occupies all sites of the crystal with equal probability in which case the density gradient vanishes and the macroscopic diffusion ceases.

Because the I-v exchanges are restricted only to NN sites, the rate of change of the impurity density on site is determined solely by the I-v configuration on this site and its NNs. Thus, the average impurity density at site may grow when the site is occupied by a vacancy while the impurity occupies one of the NN sites, so that the impurity can jump at the vacant site by the I-v exchange and the density may diminish when site is occupied by the impurity and there is a NN vacancy for the impurity to leave the site. The exchanges occur with frequency and their rates are proportional to the joint probability density of simultaneous presence of the I-v pair at the NN sites as


where is the set of 12 vectors connecting NN sites in the FCC lattice and is the vacancy occupation number. The joint pair probability density is defined as


where is the relative position of the vacancy with respect to the impurity at site ; the time dependence on the right hand side (r.h.s.) is implicit in the nonequilibrium average Eq. (97).

Thus, as could be expected on the basis of the general rate equation Eq. (98), the time derivative of the average impurity density in Eq. (3) depends on the average densities in Eq. (4) so to solve Eq. (3) further REs need be derived. Fortunately, it is easy to see that to the first order in the vacancy concentration the chain of the REs terminates already at the equation for the joint probability density . Indeed, because there is only one impurity in the system, higher order correlation functions may contain only additional vacancy occupation numbers. But then they will be of higher order in and thus can be omitted in the approximation.

The second (and the last in the approximation chosen) RE reads


It is derived as follows. By the chain rule, the time derivative on the r.h.s. of Eq. (4) should consist of two contributions: one corresponding to the change in the impurity density and the other one in the vacancy density. Therefore, the first line on the r.h.s. of Eq. (II) should coincide with the r.h.s. of Eq. (3) only corrected by the factor to account for the fact that the impurity can move only when the only available vacancy at is at a NN site to the impurity. Similarly, the second line describes the diffusion of a vacancy and formally has the same structure as the equation for the free vacancy diffusion in Appendix B but with an important difference that the jump frequencies are not all equal to but depend on the vacancy position with respect to the impurity and on the jump direction, as defined in the 5FM (see, e. g., Fig. 1 in I). In matrix notation it can be described with the use of matrix having the same structure as matrix in Eq. (101) but instead of at the matrix diagonal in should contain the sum of all jump frequencies from site to the 12 NN sites (the negative term on the last line of Eq. (II)) while the positive term should include frequencies of all possible vacancy jumps from the sites that are NN to to this site (in they are all equal to ).

As is seen, the equation for the rate of change of the pair density Eq. (II) in approximation depends only on the pair density itself. So the density evolution can be found from this equation alone, provided the initial condition is known. With the initial location of the impurity at the coordinates origin, it remains to chose the initial distribution of the vacancies. In principle, in an out-of-equilibrium system the distribution can be arbitrary. But in order to describe the evolution towards thermal equilibrium one would need to know in detail the vacancy kinetics. For example, if the vacancies are in excess in comparison with the equilibrium concentration than in order to describe such phenomena as the creation of divacancies, vacancy pores, or dislocation loops one would need to know intervacancy interactions. Besides, both under the vacancy excess or deficit their sources or sinks must be introduced into the model in order to correctly describe the kinetics leading to the equilibrium. Such terms, however, are absent in the 5FM, so the model implicitly presumes the vacancy concentration to be a constant parameter independent of the kinetics. This is possible only at thermal equilibrium. This reasoning, however, is valid only for the global vacancy concentration while a local distribution in the vicinity of the impurity still can be arbitrary. The local density perturbation concerns only number of vacancies which cannot influence global concentration in the thermodynamic limit. But to farther simplify the problem, in the present paper we will restrict our consideration to the simplest symmetric distribution. Namely, we assume that at the vacancies are distributed homogeneously everywhere in the crystal with the equilibrium density except at the site occupied by the impurity where the density should be zero by the vacancy definition and at the sites NN to the impurity where the concentration is assumed to be different from because, e. g., of the I-v interaction. These conditions are satisfied by the expression


where the first factor on the r.h.s. accounts for the initial impurity position.

The initial value problem for a linear equation with constant coefficients can be conveniently solved by means of the combined integral Laplace and Fourier transform (the LF-transform) that for an arbitrary function can be defined as


Under the transform Eq. (II) becomes




The first line on the r.h.s. in Eq. (II) is the Fourier transformed initial value of the pair density Eq. (6). The remaining two lines originate from the first and the second lines on the r.h.s. of Eq. (II). Thus, for each we obtain an infinite inhomogeneous system of equations for the -dimensional vector with components


which in the thermodynamic limit becomes infinite-dimensional. In the next section it will be shown that if the jump rates in Eq. (II) differ from the bulk rate only within a finite region around the impurity, the infinite set of equations can be reduced to a finite linear system.

To conclude this section let us derive an expression for the impurity GF in terms of the solution of Eq. (II). The LF-transformed Eq. (3) reads


where use has been made of the fact that summations over and are equivalent due to the lattice symmetry. Casting Eq. (11) in the form




is the diffusion kernel, one can recognize in Eq. (12) the first two terms of the expansion of of the GF in the Dyson form


as discussed in Sec. IV in I. In the theory of vacancy-mediated diffusion Eq. (14) can also be obtained in the rate equation approach.Bender and Schroeder (1979) In Appendix C we show that it can be rigorously derived in the framework of the Mori-Zwanzig memory-function formalism. But the equations for the memory function (or the self-energyForster (1975)) that can be obtained from Eq. (124) become useful mainly in high order expansions where they significantly reduce the number of contributions by restricting them only to irreducible ones. The lowest order term that we are interested in is easily recoverable from the much simpler REs for GF.

Iii Solution of the rate equations by the Cramer’s rule

The method of solution presented below is a generalization to the case of the 5FM of the method suggested in Refs. Tahir-Kheli and Ellott, 1982; Tahir-Kheli and Elliott, 1983 for the cases of self-diffusion and for the impurity diffusion in a two-frequency model. As will be seen, the method can be applied to any model with I-v interactions of finite range. Its essence is the reduction of the infinite (in the thermodynamic limit) set of equations Eq. (II) to a finite linear system.

To begin with, we note that one component of can be found immediately. Examination of Eq. (II) shows that satisfies the equation independently of the values of other components of . Indeed, only the terms on the third line on the r.h.s. of Eq. (II) could potentially cause problems. But because simultaneous presence of the vacancy and the impurity at the same site are forbidden, the jump rates to- and from the impurity site must be set equal to zero


Here the fictitious sixth frequency was introduced in order to formally treat the jumps between all sites in the lattice on the same grounds. Thus, satisfies Eq. (II) and later will be excluded from consideration.

Our next step is to subtract from both sides of Eq. (II) vector with the components


where is the matrix from Eq. (101) that describes the vacancy diffusion in the pure host crystal. The last line in Eq. (II) has similar formal structure and can be written down as so that after the subtraction the terms on the line become


Because beyond the forth coordination shell (CS) the rates in are equal to those of (see the discussion following Eq. (II)), the matrix difference in Eq. (17) has nonzero matrix elements only for in a vicinity of the impurity. This means that Eq. (17) contains only a finite number of components of . Taking into account that on the second line in Eq. (II) are restricted to the first CS, the number of the vector components on the r.h.s. of Eq. (II) is also finite.

As can be seen from the definition of the 5FM and from Fig. 1, there are 54 sites in four CSs plus one site at (55 in totalKoiwa and Ishioka (1983)). In matrix notation the equation can now be written as


where in the matrix we gathered all matrix elements of Eq. (17) and also the terms on the second line of Eq. (II). The inhomogeneous term is the vector composed of the terms on the first line of Eq. (II)

Figure 1: Axisymmetric diffusion along direction of the FCC lattice as described by the 5FM. The lattice sites in the vicinity of the impurity (black circle) belong to 13 equivalence classes marked as: Black points—the classified sites in the vertexes of large cubes; crosses—sites at the centers of the external faces of the cubes in the drawing, open circles—on the internal faces.

However, on the l.h.s. of Eq. (18) all components of the vector are still present after the subtraction:


The components are mutually coupled because of the structure of matrix in Eq. (101) which always connects the vector components in a given CS to higher CSs. To overcome this difficulty, let us multiply both sides of the equation by the matrix


with the matrix elements given by Eq. (106). Under the multiplication the number of unknown components of on the r.h.s. will not change but on the l.h.s. the components now decouple from each other


So by retaining only those equations that contain on the l.h.s. the same vector components that are present on the r.h.s. one arrives at a finite system of linear equations which in the conventional form reads


where is the unit matrix. The size of the system can be restricted to 54 equations because the component is already known. As is easy to see, to obtain the matrix for the reduced equation it is sufficient to multiply matrix by matrix by omitting the row and the column corresponding to in the initial matrices and , respectively. The free term on the r.h.s. can be found explicitly with the use of Eqs. (II), (106), and (107) as


With the matrix on the l.h.s. and the inhomogeneous term given by Eq. (24) being known explicitly, Eq. (23) can now be solved by the standard means of linear algebra. From Eqs. (13) and (14) it is seen that for the calculation of the impurity GF for general Fourier momentum only twelve NN components of are needed. But because FCC lattice is centrosymmetric this number can be reduced to six, so the use of the Cramer’s rule seems to be computationally simpler than the full matrix inversion.

The GF obtained from the solution of the 54 equations Eq. (23) can be used to describe situations where large values of momentum are of interest,Vogl (1996); Leitner et al. (2009); Mantl et al. (1983); Steinmetz et al. (1986) which will be illustrated in the next section.

Iv Diffusional broadening of the Mössbauer line

An important advantage of the rigorous solution for the impurity GF is that it is not restricted to small values of the Fourier momentum as the phenomenological GF but is valid for arbitrary . This makes possible theoretical description of such techniques as the quasielastic neutron scattering, the coherent X-rays and the Mössbauer spectroscopies.Vogl (1996); Leitner et al. (2009); Mantl et al. (1983); Steinmetz et al. (1986)

As an illustrative example let us consider the diffusional broadening of the Mössbauer line in the FeAl system studied in Ref. Mantl et al., 1983. The system is interesting from the point of view of the pair diffusion because of the rather strong I-v attraction of 0.29 eV as estimated from the experimental data.Mantl et al. (1983) However, the diffusion profiles simulated on the basis of the extracted 5FM parameters were found to be roughly Gaussian (see Fig. 3 in Ref. Mantl et al., 1983). Furthermore, some other conclusions drawn by the authors do not agree with the results obtained in I, so below we will attempt to clarify these issues.

Diffusional behavior contributes to the Mössbauer broadening through the van Hove correlation function which at small impurity concentrations coincides with the single impurity GF calculated in Sec. III. In the Mössbauer studies it was found that the line broadening can be adequately described in the framework of the EMWolf (1977); Bender and Schroeder (1979); Vogl (1996); Mantl et al. (1983) so we first establish connection of our approach with the EM. To this end we note that an important parameter of the EM is the average number of steps that the impurity makes during one I-v encounter. Physically is akin to the mean diffusion distance of the phenomenological theory and establishing a formal relationship between the parameters should facilitate comparison between two approaches.

Formally the diffusional contribution to the Mössbauer line broadening is given by the real part of the impurity GF at , where is the natural width of the Mössbauer line and is the energy transferred to the system by the gamma ray, so our solution for the GF should be sufficient for the task.Bender and Schroeder (1979); Mantl et al. (1983) However, in Refs. Bender and Schroeder, 1979; Mantl et al., 1983 it was pointed out that the solutions of the kind of our Eq. (14) cannot be directly compared with experimental Mössbauer spectra because the measured quantity is the width of the line which can be accurately fitted by the Lorentzian distribution. But the line shape in Eq. (14) is not Lorentzian because of the -dependent diffusion kernel . The difficulty is overcome by the EM that describes the diffusion as a sequence of repeated I-v encounters.Wolf (1977); Bender and Schroeder (1979) In the course of the encounter a vacancy is always present in the vicinity of the impurity so all impurity jumps occur within a short time interval which is small in comparison with the time interval between the encounters that is equal to the association time


So to a good approximation can be neglected on the scale of and the impurity transfer during the encounter from the initial to the final position may be considered as instantaneous. This picture can be translated into the following approximation to GF. The diffusion equation satisfied by GF in Eqs. (14) and (123) under the inverse LF-transform would read


Because varies on the time scale which is much shorter than , the self-energy will differ from zero only in a narrow region of its argument where so in the argument of in Eq. (26) can be neglected: . Now the integration over reduces to multiplication of GF at time by where the upper limit can be safely extended to infinity.Bender and Schroeder (1979) But such integration is equivalent to component of the Laplace transform of , so the LF-transformed GF in this approximation will be


(the correction due to the natural line width was found to be negligible in the experimental conditions of Ref. Mantl et al., 1983). In this approximation the real part of GF becomes a Lorentzian of widthBender and Schroeder (1979)


and is a real non-negative function of .

The quantityMantl et al. (1983)


with corresponding to the Mössbauer gamma rays and with from the Cramer’s rule solution of Sec. III has been calculated and smeared with the experimental resolution. In Fig. 2 the curve thus obtained is compared with the experimental data and the Monte Carlo (MC) simulations of Ref. Mantl et al., 1983. The same parameters of the 5FM were used in the calculations with a minor correction for which was reduced by 4% to account for the fact that in the MC simulations the correlation factor was estimated as while Eq. (56) gives . So was modified for the calculated diffusion constant agreed with the fit to experimental data.Mantl et al. (1983)

Figure 2: Symbols and the dashed line: experimental data and the calculated anisotropy of the Mössbauer resonance from Ref. Mantl et al., 1983; solid line: Eq. (29) with the same parameters as in Ref. Mantl et al., 1983 except the value of being reduced on 4%, as explained in the text.

In the EM the momentum-dependent line broadening was found to be (see Eqs. (6) and (7) in Ref. Wolf, 1977)


where is the Fourier transform of the distribution of the impurity density over the lattice sites after the I-v encounter with being the average time between the encounters. In the phenomenological approach similar expression can be obtained by substituting Eq. (127) from Appendix D into Eq. (28):


In the phenomenological approach is equal to the association time so comparing Eqs. (31) and (30) one may conclude that


The equality here may be achieved only for asymptotically small in the region of validity of the phenomenological theory. In the Mössbauer experiment, on the other hand, is finite and comparatively large so in Eq. (32) cannot be found in the small- limit. However, below we will see that the number of impurity jumps depends only on the small momenta, Eq. (32) should be sufficient for establishing an exact relation between and .

The relationship can be found with the use of the expression


which is a slightly rearranged Eq. (8) from Ref. Wolf, 1977 written in our notation. In the thermodynamic limit can be found from the continuum inverse Fourier transform as


where dimensionless momentum was introduced to simplify notation and the use has been made of the cubic symmetry. Substituting this into Eq. (33) and using the identity


it can be seen that depends only on the behavior of at small . Indeed, the sum over in Eq. (33) is calculated by first applying the Laplacian to Eq. (35) and then carrying out the integration over by parts twice. The second derivatives at are calculated with the use of Eq. (32) to give


where we replaced general by from Eq. (131) because the phenomenological approach was developed for the case of strong I-v binding. Now substituting from Eq. (128) into Eq. (36) and using the definition of Eq. (129) one finds


where use has been made of the definition of Eq. (130). As will be discussed below, this expression does not agree with the qualitative conclusions about the dependence of on the 5FM frequencies reached in Ref. Mantl et al., 1983 on the basis of the MC simulations. The comparison may be not warranted because Eq. (37) has been derived for the strong-coupling case while simulations covered also other cases.

To clarify this issue let us calculate for the general 5FM within the approach of Ref. Bender and Schroeder, 1979. In this approach the mean number of the impurity jumps is calculated as


where is the probability for the vacancy to return from the first CS on the impurity site. can be found as follows. The rate of the vacancy return on the impurity site is equal to the rate of I-v exchanges irrespective of the binding. The rate of the definite departure of the vacancy away from the impurity is equal to the product of the total rate of the vacancy jumps from the first CS to higher shells multiplied by the probability to diffuse infinity far from the impurity . Each factor in the product can also be calculated for any values of the frequencies, not only for the strong binding case. With all rates being known, the return probability is found as the ratio of the return rate to the total rateBortz et al. (1975)


Substituting this into Eq. (38) one finds


This expression differs from Eq. (37) only on a numerical constant equal to unity while the dependence on the frequencies is the same. It is interesting to note that being derived in a more general case Eq. (40) looks as more reliable than Eq. (37) restricted to the case of tight binding. However, the latter seems to be more physical because when (hence, ) vanish it predicts the correct number of impurity jumps equal to zero while Eq. (40) predicts one jump. Nevertheless, below we will use from Eq. (40) for consistency with Ref. Mantl et al., 1983.

From Fig. 3 it is seen that from Eq. (40) agrees with the MC simulations for all sets of the 5FM parameters studied in Ref. Mantl et al., 1983. But the qualitative description of the frequency dependence of suggested in the paper does not agree with Eq. (40). First disagreement concerns the conclusion that strongly depends on ratio and weakly depends on the ratio .Mantl et al. (1983) But Eq. (37) does not depend on at all while through it strongly depends on ratio, especially when it is large, as was shown with the use of MC simulations in I and will be further confirmed below within the rigorous solution. Second, when and so in Eq. (37) will also strongly depend on the ratio while in Ref. Mantl et al., 1983 it was concluded that does not depend on it. One source of the discrepancies seems to be due to the fact that the case of large ratios was not investigated in Ref. Mantl et al., 1983 while in our study it was of major importance because the first-principles calculations discussed in I predict that it dominates strong I-v binding, at least in the aluminum host.

Figure 3: The number of impurity jumps during single I-v encounter as obtained in MC simulations in Ref. Mantl et al., 1983 and as calculated with Eq. (40); the dashed line corresponds to strict equality between the two values.

Another apparent disagreement concerns roughly Gaussian diffusion profiles obtained in the MC simulations in Ref. Mantl et al., 1983 despite the large I-v attraction which according to our findings might lead to NGDPs. This, however, is explainable by the small value of the parameter that can be calculated with the use of Eqs. (36) and (131) on the basis of the experimentally fitted frequencies and .Mantl et al. (1983) The small values of both quantities are mainly due to the small value of found in the fit. Thus even in the case of strong binding the mobility of impurity can be low because of the reduced frequency of the I-v exchanges that define the diffusion. characterizes the extent of the exponential part of the NGDP. But with being smaller than even the NN distance non-Gaussian behavior is hardly detectable in the 3D case, though in 2D the microscopic profiles were reliably measured in Refs. van Gastel et al., 2011; Grant et al., 2001. The small value of and in the FeAl system were partly due to the high experimental temperature of 923 K. There is not enough data to decide whether NGDPs can be observed in this system at lower temperature, though usually grows when the temperature is lowered.Cowern et al. (1990); Duffy et al. (2003); Tokar (2018)

V Diffusion along a symmetry axis

A major goal of the present study is to give a rigorous justification to the phenomenological theory of impurity diffusion in the FCC host developed in I. The theory has been based on the notion of the mobile intermediate state of the impurity introduced in Ref. Cowern et al., 1990 which in the case of the 5FM has been identified in I with the tightly bound I-v pairs. This has made possible derivation of simple formulas for the impurity GF and the diffusion profiles in stark contrast with the complicated expression for the rigorous solution of Sec. III. So our aim in this section is to show that the phenomenological expressions do agree with the rigorous approach in the limits of strong I-v binding and the macroscopic diffusion that corresponds to small values of and in the frequency-momentum space or, equivalently, to macroscopic spatiotemporal scale.

In cubic crystals diffusion is isotropic in the limit of small so instead of a general momentum that requires solution of 54 equation, the momentum along a symmetry axis can be chosen in order to reduce the system size. To this end let us chose the direction along axis with . In this case the components of vector in Eq. (23) can be divided into 13 equivalence classes, as shown in Fig. 1. All components within a class are equal so it is sufficient to retain only one equation for each class which reduces the system of 54 equations to a system of 13 equations


where on the r.h.s. only 13 components of vector Eq. (24) corresponding to different classes has been kept, is the identity matrix of size 13 and and are and matrices, respectively, that are obtained from and as follows. Because all vector components belonging to one class are the same, their precise coordinates are irrelevant and can be characterized only by the class number which is used as the subscripts of matrices and . To find one can take any site from class and then sum over all from class . For example,


as can be seen from Fig. 1 by taking the leftmost site from class 5 and checking that there is exactly one NN site at belonging to class 1 (the leftmost in the figure from this class), one at (the rightmost in the class) and two sites at . It is to be noted that matrices and are not symmetric. Besides, one would need also the matrix elements , where 0 refers to the “class” consisting of the site at the impurity position. are necessary to fill the 0-th column of matrix . All expressions of the matrix elements of in terms of needed in calculations below can be found in ancillary file in Ref. SM, .

Matrix can be obtained from in a similar way by summing all contributions into matrix elements between the classes:




Matrix has the familiar gain-loss structure (see Appendix A) with the barred frequencies in the off-diagonal matrix elements having plus sign and the diagonal terms with the minus sign, but not necessary positive or negative values because of the subtraction of . Thus, for example, is equal to because there is exactly four sites that are NN to a site in 4th class from which a vacancy may jump at this site with the rates different from and similarly is equal to because there is four sites (all from class 1) where the vacancy can jump from a site in class 4 (see Fig. 1).

Because our goal is to find the impurity GF, we are interested only in those components of vector that enter into Eq. (11), that is, only in the components belonging to the first and the third classes. Moreover, because for any NN vector from the first class there exists an NN vector in the third class such that


it is easily seen that one can express the solution for in terms of only one of the components because the FCC lattice is centrosymmetric. Because the initial value Eq. (6) was also assumed to be centrosymmetric, the kinetics will preserve the symmetry, so during the whole system evolution the following equality will hold


Applying the LF-transform Eq. (7) to this equation and changing the summation on the r.h.s. from to one arrives at the equality


Thus, in view of Eq. (46), only two functions with being either 1 or 3 are needed for the calculation of the impurity GF in Eq. (11) . Choosing with the use of the Cramer’s rule one gets




is the determinant of the system in Eq. (41) and is this determinant with the first column replaced by the r.h.s. vector Eq. (74) according to the Cramer’s rule. Due to the axial symmetry Eq. (13) simplifies to


where the summation is over the upper and the lower signs in the summand and is any of the four lattice vectors from the first class. Eq. (51) together with Eq. (14) solves the problem of finding the impurity GF in the uniaxial geometry.

Vi Diffusion of tightly bound I-v pairs

In the axisymmetric case the impurity GF is fully determined by the diffusion kernel Eq. (51) which in the limit in FCC lattice can also describe macroscopic diffusion for any with small absolute value.

In the phenomenological theory only the macroscopic diffusion can be treated. The corresponding diffusion kernel can be found from Eqs. (44) and (53) in I as


where was replaced by to facilitate comparison with the above rigorous expression. Thus, in order to prove equivalence of the phenomenological approach and the rigorous theory in the diffusion limit in the case of strong I-v binding it is necessary to show that Eq. (52) approximately reproduces Eq. (51) at small values of , , and (see Eq. (126)).

To this end let us first compare the values of the diffusion constant in both approaches. Using Eq. (14) it is easy to show that the diffusion constant is obtained from the diffusion kernel as


Applying this to Eq. (52) one gets


where the diffusion constant was supplied by the subscript to remind that the phenomenological expression is valid only in the case of strong I-v binding. Eq. (54) is satisfied in the phenomenological theory (see Eq. (55) in I) and we are going to show that this is also the case in the rigorous approach. This amounts to showing that agrees with the known expressions for in the tight-binding limit. But this would also hold if is correct for all values of frequencies, not only in the limiting case. Though rigorous analytic expression for in the canonical version of the 5FM has been derived in Ref. Koiwa and Ishioka, 1983, below we present its numerical calculation which has an advantage of admitting generalization to the extended versions of the model.Claire (1978); Wu et al. (2016); Bocquet (2014); Tuijn et al. (1992)

vi.1 The diffusion constant

Correlation effects are conventionally accounted for in the diffusion constant through the correlation factor defined asClaire (1978); Manning (1968); Philibert (1991)


where in the 5FM case was shown to have the general formManning (1968); Koiwa and Ishioka (1983); Bocquet (2014)


This expression has been thoroughly studied in literature and so may serve as a stringent test for new calculation techniques.

The limit in Eq. (53) can be calculated by first expanding the -dependent quantities in the nominator and denominator of Eq. (51) to the second order in as


In this notation


The limit is straightforward to calculate for the first term on the r.h.s. of Eq. (57). As is easily seen from the definition of the LF transform Eq. (7), if the correlation function Eq. (4) has a non-zero asymptotic when , the Laplace part of the transform develops the singularity . But at large times and in the absence of macroscopic fluxes, which is the case when , all kinetics die out and the correlation function tends to its equilibrium value


where is the equilibrium vacancy concentration at site relative to the impurity; the factor appears due to the fact that as the only impurity in the system occupies any of sites with equal probability . The LF transform Eq. (7) of Eq. (59) at gives


Thus, the first term on the r.h.s. of Eq. (57) in the limit is equal to


where use have been made of Eqs. (125) and (126).

The limit in the last term of Eq. (57) concerns only so it is convenient to first obtain it with the use of Eq. (24)


and (49) as


where differs from in that instead of the limiting value Eq. (62) is substituted and, besides, the vacancy concentration is factored out of the determinant so that the first column of should be replaced simply by the column of unities. Further, using the reflection symmetry of the 1D geometry it can be shown that determinant is an even function of , so only contributes to the derivative in Eq. (57). The dependence of on comes only from the matrix element


of matrix Eq. (V) that contributes to the third column of matrix the terms of the form:


Similar contributions to column one from the matrix element disappear from because this column is replaced by unities. Thus, differentiation of with respect to in the limit amounts to the following expression for the second term in the diffusion constant Eq. (57)


where is the determinant of matrix where the first column is replaced with unities and the third one with the terms . Now adding the two terms in Eq. (57) in the limit with the use of Eqs. (60), (126), and (66) and using the definition of the correlation factor Eq. (55) one arrives at


To compare with in Eq. (56) the latter was solved with respect to as


and from Eq. (67) was calculated for more than randomly generated quartets of the ratios , . The calculated values were substituted in Eq. (68) and as can be seen from Fig. 4, where for brevity only about 10% of the simulation data are presented, no noticeable dependence on , is discernible. This supports the results of Refs. Manning, 1968; Koiwa and Ishioka, 1983 and confirms the assumption that the general form of in Eq. (56) is exact to the leading order in .

Figure 4: Upper curve: points for from Eq. (68) calculated for random values of ratios , ; on the scale of the figure the data would be indistinguishable from the curves of Refs. Manning, 1968; Koiwa and Ishioka, 1983. Lower curve: points calculated with the use of Eq. (71) for random ratios , , ; symbols: MC simulations from I.

vi.2 Diffusion kernel in the strong-binding limit

The expression for obtained from Eqs. (55) and (67) is valid for all values of the frequencies, including the tight-binding case. Thus, Eq. (54) provides one relation between three phenomenological parameters and the known diffusion constant so in order to calculate all of them in the rigorous approach two additional relations are needed. Let us begin with the calculation of the decay rate . It does not enter as a parameter in the rigorous solution but has to be found from the assumption that in the tight-binding case the exact diffusion kernel has the form of Eq. (52) and that depends only on and as in the phenomenological expression Eq. (130). Confirmation of these assumptions in the framework of the rigorous theory would provide a nontrivial check of the phenomenological approach.

As is seen from Eq. (126), strong I-v binding means either small value of , or large value of , or both. The difficulty in studying the tight-binding limit within the rigorous approach is that the cases and should be studied separately because the frequencies enter in a different manner in both the decay rate Eq. (130) and in the determinants that define the diffusion kernel, as can be seen from Eqs. (V), (49), (41), and (13). Therefore, let us start from the case. We first note that the determinants in Eq. (52) are polynomial in their matrix elements so the singular behavior at may originate only from the determinant in the denominator of Eq. (49) which in this limit should be representable as


where and are underlined to indicate that in this expressions they in general are functions of and . When the latter variables are equal to zero the underlined quantities would coincide with the constants and . The function is fixed by the condition that in the second factor has the coefficient equal to unity when , , and are equal to zero and ; for brevity, only the variables that change their values in the expressions below are explicitly shown as the arguments of with other parameters being treated as constants.

From Eq. (69) one can derive the derivatives of the determinant over its arguments as


(additionally assuming that ) so that if the decay rate has the form of Eq. (130) then


should be the function of only the ratio and should agree with the runaway probability defined in I. In Fig. 4 it can be seen that this is indeed the case.



calculated for randomly generated frequency ratios agreed with Eq. (129) to within the accuracy of the calculations. With the use of Eq. (54) and taking into account that , , and in the tight-binding case have the same values as in the phenomenological theory to the accuracy of the numerical calculations, we conclude that the value of the association rate in the rigorous approach also has the correct value.

To finalize comparison with the phenomenological theory in the case it remained to consider the contribution due to the associated pairs that may be present in the initial state. The contribution is given by the term in Eq. (52) proportional to


In Eq. (51) this term should originate from Eq. (49) through the similar contribution in the initial condition Eq. (24)


In the Cramer’s rule solution it can be separated from the rest of the terms because determinants are linear functions of the column vectors. Explicit expression for can be derived from Eqs. (51) and (49) with in the last equation replaced by , where is the vector of size 13 composed of from Eq. (74) with belonging to the respective classes (all within the class are the same). Now taking into account representation Eq. (69) for , from Eq. (51) can be cast in the form similar to Eq. (73) as


where notation is the same as in Eq. (69) and in are gathered all factors not present in the rest of the equality. Next, from Eq. (73) it is easy to see that


Similar transformation of Eq. (75) leads to the following condition of the agreement between the rigorous approach and the phenomenological theory:


The limit in the expression for in terms of the determinant ratio can be found similar to the case of the diffusion constant in Sec. VI.1. But a simpler route is to calculate it numerically. For trios of the ratios , , and has been generated and Eq. (77) has been found to be satisfied with very high accuracy.

In the second case of tight-binding all calculations were performed along the above lines with a few technical modifications. From the computational point of view the main difference between the two cases is that the determinants in Eq. (49) are polynomials of tenth order in , as can be seen from Eqs. (V) and (50) so when their ratio may be hard to compute with sufficient precision. The difficulty was resolved by dividing the last ten columns in the determinants by so both the matrix elements and the determinants became bounded as and thus easily manageable numerically. The condition in the equations above was replaced by the condition and the derivative over replaced by the derivative with respect to . The calculations performed with the use of this trick confirmed asymptotic equivalence of the rigorous solution to the phenomenological theory also in this case.

Vii I-v interaction of arbitrary strength

The phenomenological theory developed in I has been based on the notion of the mobile impurity state which in the case of the 5FM has been identified with the tightly bound I-v pairs.Cowern et al. (1990); van Gastel et al. (2001); Tokar (2018) In previous sections we succeeded in showing that the expressions derived in I can be rigorously justified within the rigorous approach but only in the limit of tight binding. But this regime is obviously not universal and in different systems I-v interaction will exhibit either only a weak or no attractionSteinmetz et al. (1986) or the repulsion.Wu et al. (2016) The rigorous solution in principle covers all possible cases but in view of the physical transparency and simplicity of the phenomenological theory it would be desirable to be able to asses its applicability in the cases when the tight-binding condition is not obviously satisfied.Mantl et al. (1983) Besides, the pair diffusion picture predicts such interesting and unusual phenomena as the NGDPs and the non-Fickian diffusion but as was noted in the Introduction, in Refs. Brummelhuis and Hilhorst, 1988; Toroczkai, 1997; van Gastel et al., 2011 it was shown that NGDPs can occur even in the absence of I-v attraction. This poses the question of whether the phenomenological theory does not miss possible additional contributions into the NGDP phenomenon. This question should also be addressed within the rigorous approach.

Unfortunately, the rigorous solution is difficult to deal with because of its complexity. Even in the simpler uniaxial case it has been expressed through complex determinants which apparently can be calculated only numerically. Besides, to find the diffusion profiles the solution has to be inversely LF-transformed back to the space-time variables which is a nontrivial task even in the simple 1D geometry. Explicit expression for the transformed impurity GF reads:


where is the number of planes along the chosen direction. The contour in the integration over in Eq. (78) is defined as the vertical line passing to the right of all singularities of the integrand.Weisstein () But it can be shown that in the 5FM all perturbations relax to the equilibrium state corresponding to the zero-eigenvalue eigenmode and all other eigenvalues are negative real numbers.van Kampen (1981); Tokar (1996, 2005) Therefore, the integration contour can be deformed as shown in Fig. 5. Now using the Green’s functions property