Nonlocal diffusion currents at the nanoscale

TL;DR

This paper derives a nonlocal diffusion equation for impurity currents at the nanoscale, revealing short-time nonlocal effects and proposing experimental tests, with implications for nanostructured materials.

Contribution

It introduces a generalized nonlocal diffusion model that accounts for boundaries and variable defect concentrations, extending existing pair mechanism theories.

Findings

Nonlocal effects are significant at short diffusion times.

Predictions align with local theories at large times.

Experimental verification of uphill diffusion and impurity segregation profiles is feasible.

Abstract

An integro-differential expression for the diffusion current of the impurities diffusing by the mechanism of bound impurity-defect pairs has been derived. The ensuing nonlocal diffusion equation generalizes the existing theories of diffusion by the pair mechanism in unbounded systems with homogeneous defect distributions on the systems with boundaries and variable defect concentration. It has been established that the nonlocality manifests itself only at short diffusion times while at large times, in particular, under the stationary conditions, the predictions of the nonlocal theory would coincide with existing approaches based on local diffusion currents. Possibilities of experimental verification of the nonlocal theory have been suggested. In particular, the explanation of the uphill diffusion in silicon by the pair drag can be tested within slightly modified conventional experimental setups. Also, it has been shown that in the nonlocal theory the impurity segregation profile caused by the drag should differ at the initial stage of evolution from the predictions of the local theories. Because the segregation influences mechanical, corrosive, and electronic properties of materials, the nonlocal character of the pair diffusion may have important implications for nanostructured materials.

Equations

$$C_p(\mathbf{r},t)\propto C_i(\mathbf{r},t)C_d(\mathbf{r},t)$$

$${\mathbf{J}}_i^{(loc)}\propto -D_p\nabla C_p=-\frac{D_i}{C_d^{\ast }}{C}_i{C}_d(\nabla \ln C_i+\nabla C_d)$$

$${\mathbf{J}}_i^{(loc)}=-L_{ii}\frac{\nabla {\mu }_i}{k_{B}T}-L_{iv}\frac{\nabla {\mu }_v}{k_{B}T}.$$

$$\frac{\partial }{\partial t}{G}_p({\mathbf{r},\mathbf{r}}_0,t)=D_p{\nabla }^2{G}_p({\mathbf{r},\mathbf{r}}_0,t)-ϵG_p({\mathbf{r},\mathbf{r}}_0,t),$$

$$G_p({\mathbf{r},\mathbf{r}}_0,t)=\delta ({\mathbf{r}-\mathbf{r}}_0).$$

$$P({\mathbf{r},\mathbf{r}}_0)=ϵ{\int }_0^{\infty }{G}_p({\mathbf{r},\mathbf{r}}_0,t)dt$$

$$P({\mathbf{r},\mathbf{r}}_0){|}_{|{\mathbf{r}-\mathbf{r}}_0|\to \infty }\propto \exp (-|{\mathbf{r}-\mathbf{r}}_0|/\lambda ),$$

$$\lambda =\sqrt{D_p/ϵ}$$

$$g(\mathbf{r},t)=g^{\ast }{C}_d(\mathbf{r},t)/C_d^{\ast }$$

$$t=O(t_p)\gg \Delta t\gg \tau .$$

$$\frac{\partial }{\partial t}{C}_i(\mathbf{r},t)$$

$$-{\lambda }^2{\nabla }_{\mathbf{r}}^{2}P({\mathbf{r},\mathbf{r}}_0)+P({\mathbf{r},\mathbf{r}}_0)=\delta ({\mathbf{r}-\mathbf{r}}_0).$$

$$\mathbf{n}\cdot \nabla P({\mathbf{r},\mathbf{r}}_0){|}_B=0,$$

$$\frac{\partial }{\partial t}{C}_i(\mathbf{r},t)={\lambda }^2{\nabla }_{\mathbf{r}}^2\int P({\mathbf{r},\mathbf{r}}_0)C_i({\mathbf{r}}_0,t)g({\mathbf{r}}_0,t)d{\mathbf{r}}_0=-{\nabla }_{\mathbf{r}}\cdot {\mathbf{J}}_i(\mathbf{r}),$$

$${\mathbf{J}}_i(\mathbf{r})=-{\lambda }^2{\nabla }_{\mathbf{r}}\int P({\mathbf{r},\mathbf{r}}_0)C_i({\mathbf{r}}_0,t)g({\mathbf{r}}_0,t)d{\mathbf{r}}_0.$$

$$P({\mathbf{r},\mathbf{r}}_0){|}_{\lambda \to 0}\to \delta ({\mathbf{r}-\mathbf{r}}_0).$$

$${\mathbf{J}}_i^{(loc)}(\mathbf{r},t)\simeq -{\lambda }^2\nabla [C_i(\mathbf{r},t)g(\mathbf{r},t)]$$

$$P(\mathbf{k})=(1+{\lambda }^2{\mathbf{k}}^2{)}^{-1},$$

$$C_i(\mathbf{k},t)=\exp \left(-\frac{g{\lambda }^2{\mathbf{k}}^2}{1+{\lambda }^2{\mathbf{k}}^2}t\right)C_i(\mathbf{k},t=0).$$

$$C_i(\mathbf{r},t\to \infty )\simeq \frac{A}{g(\mathbf{r})}$$

$$\frac{\partial C_i(\mathbf{r},t)}{\partial t}={\lambda }^2{\nabla }^2[g(\mathbf{r},t)C_i(\mathbf{r},t)]\simeq \frac{D_i}{C_d^{\ast }}{\nabla }^2[C_d(\mathbf{r},t)C_i(\mathbf{r},t)]$$

$$D_i=g^{\ast }{\lambda }^2.$$

$$g(x)=\left\{\begin{array}[]{lc}g_{s}[s^{-1}+(1-s^{-1})x/d_{s}],&0\leq x<d_{s},\\ g_{s},&d_{s}\leq x\leq L,\end{array}\right.$$

$$C_i(\mathbf{r},0)=C_{i0}=Const$$

$$\bar{g}\simeq g_s\frac{\lambda }{d_s}.$$

$$C_i(x,t)\approx C_{i0}\left(1+\bar{g}t{e}^{-x/\lambda }\right).$$

$$C(\xi ,\theta )={\lambda }^{-1}\sum _{n=0}^{\infty }{P}_n(\theta ){\varphi }_n(\xi ,1),$$

$$\theta =gt,\xi =x/\lambda ,$$

$$P_n(\theta )=({\theta }^n/n!)e^{-\theta },$$

$${\varphi }_{n=0}(\xi ,1)=\delta (\xi ),$$

Full text

I Introduction

In crystalline solids the substitutional solute atoms diffuse predominantly by indirect mechanisms with the mobility of the impurities (the solute atoms in dilute alloys) induced by the mobile point defects,—the vacancies (v) and the interstitials (I).Philibert (1991); Fahey, Griffin, and Plummer (1989) Under such mechanisms the diffusion can be strongly influenced by the impurity-defect (i-d) interaction. For example, in the case of strong i-d attraction the tightly bound complexes (or pairs) formLidiard (1955); Krivoglaz (1961) which may considerably enhance the impurity mobility because the diffusion-mediating defect is always available within the pair. Besides, the pair diffusion may underlie the solute segregation at the grain boundaries in quenchedAust et al. (1968); Anthony and Hanneman (1968) and/or irradiated metals.Johnson and Lam (1976); Wagner et al. (1983); Ardell and Bellon (2016) The excess defects flowing to the boundary sinks may drag along the bound impurities and if the defect gradients are sufficiently large the impurities will diffuse uphill, that is, from small to large solute concentrations, thus causing their pileup near the boundary.

Segregation at the grain boundaries and the external surfaces is of great practical interest because it may significantly impact the mechanical and corrosive properties of metals (see Ref. Ardell and Bellon, 2016 and references therein) and also the electronic properties of semiconductor devices.Fahey, Griffin, and Plummer (1989); Ferri et al. (2007) In modeling the segregation, various theoretical techniques have been used such as the phenomenological rate equations,Johnson and Lam (1976); Wagner et al. (1983); Ardell and Bellon (2016); Fahey, Griffin, and Plummer (1989); Cowern et al. (1990), generalized Fick’s laws and the Onsager formalism,Garnier et al. (2013); Ardell and Bellon (2016) as well as the Monte CarloBarashev (2005); Soisson (2006) and the molecular dynamic simulations.Domain and Becquart (2005); Azeem et al. (2019)

The theoretical studies revealed that in the systems with high impurity concentration and highly supersaturated defects of both kinds the pair diffusion represents only one of many diffusion routes which hinders its detailed investigation in such systems. To facilitate the task, it seems reasonable to resort to model systems where (i) the pair mechanism dominates the diffusion, (ii) is driven by a small concentration of defects of only one kind: the vacancy or the interstitial, and (iii) in dilute alloys where the i-i interactions can be neglected. This idealized picture was approximately realized in Refs. Cowern et al., 1990, 1991, 1992 in experiments on diffusion of boron impurities in silicon. The simplicity of experimental setup allowed the authors to describe the diffusion with only two rate equations which made possible a detailed analysis of the pair mechanism. The most important outcome of the analysis of Ref. Cowern et al., 1990 was that the pair diffusion is intrinsically non-Fickian (n-F). The reason is that when the pair nucleates at some point in the crystal and then diffuses on a macroscopic distance away from it, the impurity current carried by the pair is defined to a large extent by the densities of the impurities and the defects at the initial point but is not related to the local concentrations at the actual position of the pair in contradiction to the first Fick’s law. The nonlocality of the pair currents was not noted in the segregation studiesJohnson and Lam (1976); Soisson (2006); Garnier et al. (2013); Ardell and Bellon (2016) though it may have important experimental consequences, for example, produce n-F exponential diffusion profiles.Cowern et al. (1990)

The aim of the present paper is to derive a phenomenological expression for the nonlocal current arising in the diffusion by the pair mechanism and to use it to extend the description of the pair diffusion developed in Refs. Cowern et al., 1990; van Gastel et al., 2001 for systems with constant defect concentration to the case of bounded systems with inhomogeneous defect distribution. This extension will make possible the description of the uphill diffusion in siliconWang et al. (2001); Ferri et al. (2007) by the same pair drag mechanism as in metals.Aust et al. (1968); Anthony and Hanneman (1968) Furthermore, it will predict the possibility of observing the exponential diffusion profiles seen in the homogeneous systemsCowern et al. (1990); van Gastel et al. (2001) in the impurity segregation profiles caused by the drag in strongly varying defect distributions.

II Diffusion by the mechanism of mobile i-d pairs

To begin with let us qualitatively discuss the description of the pair diffusion based on the local currents. In the continuum approximation the pair density at point ${\bf r}$ at time $t$ should be proportional to the local concentrations of the impurities and of the defectsJohnson and Lam (1976); Morehead and Lever (1986); Mulvaney and Richardson (1987)

[TABLE]

where the proportionality constant is given by i-d pairing rate.Morehead and Lever (1986); Mulvaney and Richardson (1987) In the continuum approximation the diffusivity of a tightly bound pair can be characterized by the diffusion constant $D_{p}$ which in microscopic models of the vacancy-mediated diffusion can be calculated explicitly.Lidiard (1955); Johnson and Lam (1976); Krivoglaz and Repetskiy (1971); Garnier et al. (2013); Tokar (2018) Because we assume that the impurity diffuses only within the pairs, the impurity current should be proportional to the pair currentMulvaney and Richardson (1987)

[TABLE]

where $D_{i}$ and $C_{d}^{*}$ are the equilibrium impurity diffusion constant and the defect concentration, respectively. As is seen, when $C_{d}({\bf r},t)=C_{d}^{*}$ the standard first Fick’s law is recovered. But when $C_{d}$ depends on the space coordinates and its logarithmic gradient is opposite to the logarithmic gradient of $C_{i}$ and, besides, exceeds it in the absolute value—the current in Eq. (2) flows along the impurity gradient causing the uphill diffusion.

In the dilute limit Eq. (2) agrees with the strong binding limit of the expression derived in the local approach in Ref. Garnier et al., 2013 for the vacancy-mediated diffusion ($d=v$). The generalized Fick’s law was expressed in terms of the Onsager matrix $\hat{L}$ asGarnier et al. (2013)

[TABLE]

Taking into account that in the dilute limit $\nabla\ln C_{i(v)}={\nabla\mu_{i(v)}}/{k_{B}T}$Yoshida (1976) and substituting this into (2) one obtains Eq. (3) with $L_{ii}=L_{iv}=D_{i}C_{i}C_{v}/C_{v}^{*}$ so that $L_{iv}/L_{ii}=1$ in accordance with the strong-binding limit in Ref. Garnier et al., 2013. It is to be noted that besides Refs. Mulvaney and Richardson, 1987; Garnier et al., 2013 the local expression Eq. (2) can be obtained from the pair contribution to the diffusion current derived in Refs. Johnson and Lam, 1976; Wagner et al., 1983; Soisson, 2006.

II.1 Nonlocal currents

A distinct consequence of the n-F nature of diffusion by the pair mechanism is that the impurity profiles observed in Ref. Cowern et al., 1990 and later in Ref. van Gastel et al., 2001 were not Gaussian as in the Fickian diffusion but had the exponential shapes. In the latter reference it was shown that this was a direct consequence of the instability of the pairs which can be shown as follows. The kernel $G_{p}$ describing the diffusion of an unstable pair decaying with the rate $\epsilon$ satisfies the equationvan Gastel et al. (2001)

[TABLE]

where ${\bf r}_{0}$ is the point at which the pair was nucleated at time $t=0$ so the initial condition for $G_{p}$ reads

[TABLE]

In a homogeneous system the solution of Eq. (4) is given by the conventional Gaussian kernel multiplied by $\exp(-\epsilon t)$ to account for the decay. As $t\to\infty$ the pair completely dissociates and the impurity atom immobilizes in the substitutional positions according to the probability distribution given by the integralvan Gastel et al. (2001)

[TABLE]

which in all physical dimensions and geometries has the exponential tailCowern et al. (1990); van Gastel et al. (2001); Tokar (2018)

[TABLE]

where

[TABLE]

is the average migration distance of the impurity during one i-d encounter.Cowern et al. (1990); van Gastel et al. (2001)

Because both the pair diffusion and the decay are determined by the movements of the defect, the kinetics described by (4) proceed at the microscopic time scale. The movements of the impurity at this sale can be so fast that sometimes even cannot be detected experimentally.van Gastel et al. (2000) In practical applications, however, one is usually interested not in the diffusion of individual impurities and/or pairs but in the evolution of the concentration profile $C_{i}({\bf r},t)$. Its evolution takes place at a slow macroscopic time scale because the speed of the profile change is limited by the concentration of lattice defects $C_{d}$ which is usually small due to their high creation energy in the range of several eV.Fahey, Griffin, and Plummer (1989); Philibert (1991) Therefore, the substitutional impurity remains immobile during long intervals between the i-d encounters which trigger short bursts of the diffusion steps.van Gastel et al. (2000) The average time span between the encounters $t_{p}=g^{-1}$, where $g\ll 1$ is the i-d pairing rate which in the local approximation is proportional to $C_{d}$Cowern et al. (1991)

[TABLE]

where the starred quantities are the equilibrium values. In the absence of i-d attraction the encounter may amount to only one impurity displacement on the distance of order of the lattice constant $a$. However, under strong attraction when the decay rate $\epsilon$ in Eq. (8) is small the bound pair may perform a large number of the elementary diffusion steps and migrate at a distance $\lambda\gg a$.Cowern et al. (1990); van Gastel et al. (2000) In this case the effective microscopic time scale is defined by the pair lifetime $\tau=\epsilon^{-1}$ and the effective diffusion step $\lambda$.

If one is interested only in macroscopic evolution at the time scale $t=O(t_{p})\gg\tau$, the infinitesimal time step $\Delta t$ in the governing equation can be chosen in such a way that it was microscopically large but still small at the macroscopic scale:

[TABLE]

For example, in the experiments in Refs. van Gastel et al., 2000 the average time interval between the jumps was $10^{-8}$ s and during the pair lifetime the impurity made on average $\sim 10$ steps so $\tau$ was of $O(10^{-7}$ s). The mean time between the pairings, on the other hand, was $\sim 10$ s, so with the choice of $\Delta t=10^{-3}$ s both inequalities in (10) are easily satisfied. Because of the four orders of magnitude difference between $\Delta t$ and $\tau$ in this example, the creations and decays of the pair will be confined within the time step $\Delta t$ with the relative error $O(10^{-4})$. So on the macroscopic time scale the redistribution of the impurity density due to the pair diffusion looks as instantaneous. This means that within $\Delta t$ the impurity distribution after one I-d encounter can be approximated by $P$ in Eq. (6).van Gastel et al. (2001)

Thus, on the macroscopic scale the impurity diffusion looks like its disappearance at point ${\bf r}_{0}$ due to the pairing with a defect followed by its immediate redistribution in the vicinity of that point according to the probability density $P$ in (6). By analogy with Fick’s second law, the governing equation for the diffusion by this mechanism can be obtained as the condition of the impurity conservation as

[TABLE]

Here the first term on the right hand side (r.h.s.) describes the loss of impurities by the concentration profile due to the pairings and the second term describes its replenishment by the decays.

Because the pair diffusion kernel $G_{p}$ enters Eq. (11) only through $P$ from Eq. (6), the formalism can be simplified by dealing directly with $P$ without the resort to $G_{p}$. This is achieved by integrating Eq. (4) over $t$ from zero to infinity with the use of the initial condition Eq. (5) which amounts to the differential equation for $P$

[TABLE]

It should be supplemented by the boundary conditions which can be derived from the definition of $P$ in terms of $G_{p}$ Eq. (6). Because the lattice defect, by definition, cannot exist beyond the lattice, the pair cannot cross the system boundary $B$ which leads to the zero Neumann boundary condition for the pair flux through the boundary for $G_{p}$. By Eqs. (6) and (2) this translates to $P$ as

[TABLE]

where ${\bf n}$ is the vector normal to the boundary.

With the use of Eq. (12) Eq. (11) can be cast in the form of the second Fick’s law as follows. According to Eq. (12), $P$ is equal to the sum of the delta-function and of the Laplacian multiplied by $\Lambda^{2}$. Substituting it in Eq. (11) one gets

[TABLE]

where the impurity current

[TABLE]

As is seen, in contrast to the local expression (2) current in Eq. (15) is nonlocal. It depends depends on the weighted concentrations of impurities and the defects (see Eq. (9)) within a region of radius $\sim\lambda$ defined by the exponential attenuation of $P$ with distance (see Eq. (7)). Physically the region corresponds to the neighborhood of point ${\bf r}$ from which the impurity may reach the point in one migration event.

When $\lambda$ tends to zero, the current in Eq. (15) becomes local, as can be seen from (12):

[TABLE]

The migration distance,however, is not a tunable parameter but has a concrete finite value specific to the system. But because it has the dimension of length, it should be compared with the typical length scale on which the functions entering the integrand in Eq. (15) vary. Assuming that on the scale of $\lambda$ the variation of both $C_{i}$ and $g$ is weak, all functions of ${\bf r}$, $P$ in the integrand of Eqs. (11) and (14) can be replaced by the local expression (16) to give

[TABLE]

which in view of Eq. (9) coincides with the local current Eq. (2). Though in explicit calculations below Eq. (2) will always be used, notation $g({\bf r},t)$ will be kept in order to remind that Eq. (9) is only an approximation and in a rigorous treatment $g$ must be calculated with the use of the diffusion equation for the defects and with an appropriate model for i-d pairing. In particular, because both the impurity and the defect are present within the pair, a contribution similar to Eq. (11) and/or Eq. (14) should be accounted for also in the defect diffusion. In the present study it is neglected because in contrast to the impurity the defects are the most mobile entities in the system so when the impurity concentration is small the pair contribution which is proportional to $C_{i}$ is negligible. This approximation, however, may not be appropriate in concentrated alloys.Mulvaney and Richardson (1987)

Eq. (14) together with Eqs. (12), (13), and (11) are the main results of the present study. Besides the heuristic derivation presented above, they agree with the rigorous theory of Ref. Garnier et al., 2013 in the limit of slow variation of the defect and the impurity concentrations as well as with the heuristic approaches of Refs. Johnson and Lam, 1976; Mulvaney and Richardson, 1987; Soisson, 2006. The equations extend the picture of the pair diffusion developed in Refs. Cowern et al., 1990; van Gastel et al., 2001 for unbounded systems with constant defect concentration to the inhomogeneous case and systems with boundaries. In particular, the nonlocal equations exactly reproduce the evolution of 1D impurity concentration profile described in Ref. Cowern et al., 1990. In an unbounded host with constant $g$ Eqs. (12) and (14) or, equivalently, Eq. (11), can be solved with the use of the Fourier transform as follows. The solution of Eq. (12) reads

[TABLE]

where ${\bf k}$ is the Fourier momentum. With known $P({\bf k})$, the solution of Eq. (14) is easily found to be given by

[TABLE]

In Appendix A it has been shown that in 1D geometry with the initial delta-function profile $C_{i}({\bf k},t=0)=1$, where $k=k_{x}$ is the $x$-component of ${\bf k}$, the inverse Fourier transform of Eq. (19) coincides with the series solution of Ref. Cowern et al., 1990. It is important to stress that this solution, hence, the underlying physical picture, were experimentally confirmed in Refs. Cowern et al., 1990, 1991, 1992; Mirabella et al., 2013 so the nonlocal equations generalize this picture to a broader class of the pair diffusion problems.

III Impurity diffusion in stationary defect profiles

To estimate the uphill diffusion and/or the segregation predicted by the nonlocal equations let us consider the evolution of an impurity profile under a stationary inhomogeneous defect distribution. The inhomogeneity naturally arises, for example, when the concentration of defects in the crystal bulk exceeds their value at the surface, so in the absence of other sinks they flow towards it to reduce the supersaturation. In the case of sufficiently strong i-d attraction the impurity dragged by the defects will diffuse in the same direction, as follows from Eq. (2).Aust et al. (1968); Anthony and Hanneman (1968); Garnier et al. (2013) As has been pointed out previously, the drag may take place even when the impurity concentration in the bulk is lower than near the surface. It is this uphill diffusion that underlies the alloy segregation by the pair drag mechanism.Aust et al. (1968); Anthony and Hanneman (1968)

The strict stationarity of the defect distribution can be achieved only in the presence of time-independent defect sources, such as the stationary external irradiationJohnson and Lam (1976); Wagner et al. (1983); Soisson (2006); Bruno et al. (2009); Ardell and Bellon (2016) or in an externally imposed constant temperature gradient.Jaffe and Shewmon (1964) But as will be seen below, of major interest to the present study is the initial stage of segregation, so the stationarity should be a reasonable approximation even in cases of transient defect supersaturation as produced, e.g., by the quench,Aust et al. (1968); Anthony and Hanneman (1968) by the ion implantation,Fahey, Griffin, and Plummer (1989); Ferri et al. (2007) by the oxidation,Cowern et al. (1990) etc. The qualitative behaviors discussed below should also hold in the case of time-dependent supersaturation because the underlying mechanism of the pair diffusion will be operative also in this case.

The time-independence of defect distributions significantly facilitate analysis because in this case the pairing rate $g({\bf r})$ entering the diffusion equations will also be stationary. Furthermore, because diffusion equations describe relaxation processes, at large times the impurity profiles should also asymptotically approach stationary shapes. With the use of Eq. (14) it can be shown that as $t\to\infty$ the impurity distribution will acquire the form

[TABLE]

where $A$ is a normalization constant defined by the total number of impurities in the system. The validity of Eq. (20) can be seen by observing that $P({\bf r,r}_{0})$ in Eq. (12) is a symmetric function of its arguments which is a property of Green’s functions satisfying the Neumann boundary condition Eq. (13). So the integral in Eq. (14) over ${\bf r}_{0}$ is equal to the integral in Eq. (12) over ${\bf r}$, i.e., is equal to unity, which makes the r.h.s. of Eq. (14) equal to zero. The peculiarity of the profile Eq. (20) is that it does not depend on $\lambda$ which means that it will be the same for both large and small $\lambda\to 0$, that is, in the local theory. In the latter case substituting Eq. (2) into the second Fick’s law Eq. (14) one arrives at the local diffusion equation

[TABLE]

where the second equality follows from Eq. (9) and from the expressionCowern et al. (1990)

[TABLE]

As is easily seen, the stationary profile Eq. (20) satisfies also the local diffusion equation (21).

Equations similar to Eq. (21) are frequently used for the description of the pair-driven drag. In particular, in the BCCGarnier et al. (2013) and in the FCC metals in the vacancy-mediated diffusion (see equation (30) in Ref. Johnson and Lam, 1976 in the strong binding limitPhilibert (1991) $w_{3}\to 0$). In the case of interstitial defects $d=I$ the last of Eqs. (11) in Ref. Soisson, 2006 or the corresponding equation in Ref. Wagner et al., 1983 would acquire the form of Eq. (21) if the vacancy contribution were neglected. The theories in these papers were intended to deal mainly with the drag under the stationary irradiation and the stationary segregation profiles so the local should be sufficient. However, in the next section it will be shown that at small diffusion times the local form may lead to erroneous conclusions about the shape of the impurity distributions so more accurate nonlocal equations should be used in the simulations of the short-time pair diffusion.

IV Impurity drag towards the surface by the pairs

To make a semi-quantitative comparison of the segregation predicted by the nonlocal diffusion equations with experimental observations and the local theory, the impurity profiles in a symmetric film of thickness $2L$ has been simulated.Johnson and Lam (1976) Our interest was in the profiles that strongly vary at distances comparable to $\lambda$ when the local approximation can break down.

The pairing rate was modeled by a simple stationary shape

[TABLE]

where $d_{s}$ is the depth at which the pair nucleation rate saturates to a constant value $s=g_{s}/g(0)$ which will be called supersaturation and assumed to be much larger than unity.Cowern et al. (1990) For large $d_{s}$ the variation of $g(x)$ in Eq. (23) will be slow so the choice of Eq. (23) can be justified by Eq. (9) as follows. If starting from depth $d_{s}$ there exist some stationary defect sources, such as the interstitial clusters,Ferri et al. (2007) which keep a constant level of defect supersaturation for $x\geq d_{s}$ and if at the surface the defect concentration drops to much smaller equilibrium value $g(0)=g^{*}$,Johnson and Lam (1976); Ferri et al. (2007) the conventional 1D diffusion equation would predict the linear defect profile in the interval $0\leq x\leq d_{s}$ and by Eq. (9) also the pairing profile.

The uphill diffusion for a system roughly modeling the boron diffusion in siliconFerri et al. (2007) has been simulated with the use of Eqs. (11) and (23). The input parameters in the simulations have been chosen to be similar to those in Refs. Cowern et al., 1990; Pichler et al., 2004; Ferri et al., 2007. The parameters entering Eq. (23) were assessed departing from the fact that in Ref. Pichler et al., 2004 strong gradients of the defect concentration were observed below the depth $\sim 10$ nm so $d_{s}$ was set at this value and supersaturation $s=10^{2}$ has been assumed as in Ref. Cowern et al., 1990; $g_{s}=sg^{*}$ was calculated from the expression for the diffusion constant Eq. (22) assuming $\lambda\simeq 5$ nm.Ferri et al. (2007) The equilibrium diffusion constant was taken from Ref. Haddara et al., 2000. The film thickness $L=200$ nm used in the simulations was chosen to be larger than in experiments of Ref. Ferri et al., 2007 in order to exclude the influence of the second surface of the film and to justify the use of the simple expression Eq. (36) for $P$ valid for $L\gg\lambda$. This should be adequate for the simulation of impurity pileup near the surface because the experimental data in this region are very similar both for the film and for the surface of the bulk material.Ferri et al. (2007)

The results of simulations shown in Fig. 1 are similar to the experimental profiles of Ref. Ferri et al., 2007 but a slight difference can be noted. Namely, the profile corresponding to the lowest temperature has slightly positive curvature near the surface which is absent on other two curves and at the limiting curve Eq. (20). This can be attributed to the fact that with the same diffusion time (60 s) the smaller temperature means earlier stage of diffusion so the shape has not yet reached the asymptotic regime. In experiments of Ref. Ferri et al., 2007 all three curves are qualitatively similar which means that they correspond to the advanced stage of the evolution. Thus, in experimental conditions of Ref. Ferri et al., 2007 one cannot distinguish between the local and non-local cases but our order-of-magnitude estimates show that the conditions are quite close to the preasymptotic regime where this distinction could be detected. And it seems that only the nonlocality and, more importantly, the measurement of $\lambda$ and its comparison with independently obtained value could definitely prove the pair mechanism of the uphill diffusion in B-Si system.

However, the simulations also have shown that the experimental setup of Refs. Wang et al., 2001; Ferri et al., 2007 is not very convenient for the quantitative study of the pair diffusion. The cause of difficulty lies in the initial impurity profile which is peaked near the surface (see, e.g., Ref. Ferri et al., 2007) and due to inevitable experimental errors in its measurement the additional pileup would be difficult to measure with sufficient accuracy to distinguish between theoretical predictions of different models. To reliably measure the diffusion profile it seems more natural to chose the initially homogeneous impurity distribution

[TABLE]

and also to shift $d_{s}$ further from the surface in order that the large second derivative of $C_{d}$ at $d_{s}$ (see Eq. (21)) did not obscure the pair contribution.

To this end in the simulations below the boundary between the regions of the linear and the constant behaviors in Eq. (23) was shifted towards $d_{s}=100$ nm but the gradient of the supersaturation was kept the same as previously, that is at 10 nm*-1*, which means supersaturation $s=10^{3}$ at $d_{s}$. Besides, larger value of $\lambda=10$ nm has been usedCowern et al. (1990) to simplify visualization of the profiles in the region of interest (see Fig. 2).

In Fig. 2 are shown the impurity profiles near the surface simulated with the use of Eqs. (11) and (21) with the same diffusion constant $D_{i}$ and $g(x)$ in Eq. (23) in both cases. According to Ref. Cowern et al., 1990, a natural time scale for the evolution of the diffusion profile under the pair mechanism is the number of the migration events of the i-d pairs. When the number of migrations is equal to unity the diffusion amounts to a single i-d encounter with the diffusion profile acquiring characteristic exponential shape. Large number of migrations will lead to the conventional Gaussian behavior.Cowern et al. (1990) However, under the inhomogeneous defect distribution the pair nucleation rate depends on the distance to the surface so the probability of the encounter depends on the position of the impurity within the film so the calculation of the parameter $gt$ introduced in Ref. Cowern et al., 1990 becomes ambiguous. It seems reasonable to define instead a characteristic pair nucleation rate $\bar{g}$ in the near-surface region with the use of Eq. (23) as the rate at the distance $\sim\lambda$ from the surface

[TABLE]

The diffusion profile for small values of $\bar{g}t\lesssim 1$ calculated to this accuracy in Appendix B in Eq. (37) is

[TABLE]

This expression approximately agrees with the numerical solution for $\bar{g}t=1$ (see Fig. 2). Noticeable is the exponential contribution characteristic of the pair diffusion.Cowern et al. (1990); van Gastel et al. (2001)

To visualize the details of the diffusion profiles in the narrow region near the surface, the logarithmic scales have been used in Fig. 2. Because of this, the values at $x=0$ are not shown in the plots but from the bullets that mark the half width of distributions it can be seen that the concentrations calculated with the use of Eq. (21) are twice as large at $x=0$ as at $x=0.1$ nm while the profiles simulated with the pair-diffusion mechanism Eq. (11) at $x=0$ exceed $x=0.1$ nm values only on a few percent. It is to be noted that the equation has been integrated with the step $\Delta x=0.1$ nm because anyway the continuous approximation cannot be valid at distances smaller than the atomic radii of order $O(0.1~{}\mbox{ nm})$. Thus, though at large diffusion times both mechanisms will produce the same profile Eq. (20), at small to intermediate times the profiles simulated with the local diffusion equation will be much more peaked near the surface. For example, for the pair diffusion mechanism the impurity density at the surface is an order of magnitude smaller in the cases $\bar{g}t=1,5$ and about five time smaller for $\bar{g}t=20$ (at large times the distributions converge). This may have important practical consequence for the properties of the near-surface layers. For example, the onset and the distribution of the impurity precipitates may significantly differ in the two cases.

V Conclusion

In this paper the problem of impurity diffusion by the mechanism of mobile i-d pairsCowern et al. (1990) has been considered theoretically and a phenomenological integro-differential expression for the impurity diffusion current has been derived on the basis of a simple physical picture of the pair diffusion. The expression is a nonlocal generalization of the first Fick’s law. Being substituted into the second law it leads to a nonlocal diffusion equation. Its predictions differ at the early stages of the diffusion from those of the local equations that are conventionally used in the studies of the impurity segregation.Johnson and Lam (1976); Wagner et al. (1983); Soisson (2006); Garnier et al. (2013)

The difference between predictions of the two approaches has been illustrated by simulations of the uphill diffusion and the ensuing impurity drag toward the surface under the conditions of defect supersaturation in the material bulk.Aust et al. (1968); Anthony and Hanneman (1968); Wang et al. (2001); Ferri et al. (2007) It has been shown that segregation of the impurities near the surface produced by the diffusion of bound i-d pairs has some peculiarities that may have important practical consequences. The spatial scale at which the peculiarities are the most pronounced is determined by the average migration distance of the bound i-d pair $\lambda$, which experimentally was found to be in the nanometer range in semiconductor materialsCowern et al. (1990, 1992); Ferri et al. (2007) and theoretically estimated to have similar values in metals.Tokar (2018) At the macroscopic scale much exceeding $\lambda$ it does not matter whether the elementary diffusion step is of the order of the lattice constant or is equal to $\lambda$. But if the length scale of interest is similar to $\lambda$, both quantitative and qualitative differences can be detected. The underlying cause of the differences is that at the microscopic scale the pair diffusion does not obey Fick’s lawsCowern et al. (1990) which is reflected in the nonlocal nature of the diffusion current. The n-F features originate from the permanent presence of the diffusion mediating impurity within the pair which makes possible its autonomous migration at $O(\lambda)$ distances independently of the ambient defect and impurity concentrations.Cowern et al. (1990); Bruno et al. (2009) This, however, takes place only during the pair lifetime which restricts the n-F behavior to the length scale of $O(\lambda)$. The latter is temperature dependent and at low temperatures may grow on orders of magnitudeCowern et al. (1992) which can be used as an additional means of control of the segregation profiles.

Experimentally the difference between the profiles produced by the pair diffusion and by the mechanisms that admit description by local diffusion equationsJohnson and Lam (1976); Wagner et al. (1983); Soisson (2006); Garnier et al. (2013) can be observed only at early stages of evolution because at late stages both mechanisms lead to identical stationary and/or equilibrium profiles. It should be stressed that the stationary defect distributions have been used in the present paper only for simplicity and that the underlying microscopic mechanism should be operative for time-dependent distributions as well. Short-lived supersaturation can be even advantageous because it may be used to arrest segregation at an early stage to obtain a concentration profile of desired width.

The most promising experimental setup for verification of the predictions about the pair drag in semiconductors seems to be similar to that used in Ref. Ferri et al., 2007 in the study of the uphill diffusion only the initial impurity distribution near the surface should be much flatter in order to facilitate the measurement of the small (due to the short diffusion time) differences between the profile evolution predicted in different theories, such as those discussed in Ref. Ferri et al., 2007.

The stationary irradiation conditions in metallic systems Johnson and Lam (1976); Wagner et al. (1983); Soisson (2006); Ardell and Bellon (2016) are not well suited to investigation of the nonlocality effects, in particular, because of the presence of several modes of diffusion propagated by different defects. But the quench experiments similar to those of Refs. Aust et al., 1968; Anthony and Hanneman, 1968 for solutes with strong binding to vacanciesTokar (2018) could be suitable for the purpose except that the spatial extent of the segregation profiles should be reduced to the nanoscale range.Ardell and Bellon (2016) A transient nature of the supersaturation in such experiments may turn out to be even advantageous because the pair contribution is best visible at the early stages when only a few i-d migration events occurred and this number can be controlled by the strength and the duration of the non-equilibrium defect flux.

In conclusion it should be noted that the phenomenological theory of the pair diffusion developed in the present paper is only a particular case of more general notion of the mobile state suggested in Ref. Cowern et al., 1990. The mobile state of impurity should not necessary be associated with the pairing but may be caused, for example, by the kick-out mechanism.Cowern et al. (1990, 1992) Other possibilities can be envisaged and if the transition to the mobile state is driven by the defects the formalism developed in the present paper should apply because the pairs have been assumed to be unstructured unstable point particles that can be associated, for example, with the interstitial state invoked in Refs. Johnson and Lam, 1976; Wagner et al., 1983; Soisson, 2006 to describe the radiation-induced segregation. The interstitials are unstable due to recombination processes so if the value of the mean migration path $\lambda$ is of appreciable value their isolated motion at distances of order $\lambda$ can be described by the nonlocal diffusion equations derived in the present paper. In this respect the results of the present paper may be viewed as an alternative approach to the description of the diffusion via a fast intermediate suggested in Refs. Cowern et al., 1990; Chen et al., 2010.

Acknowledgements.

I am grateful to René Monnier for carefully reading a preliminary version of the manuscript and for useful suggestions on its improvement. Also, I would like to express my gratitude to Hugues Dreyssé for support and encouragement.

Appendix A Comparison with the diffusion profile of Ref. Cowern et al., 1990

In slightly modified notation, the evolution of 1D impurity profile starting from the delta-function distribution at $t=0$ can be described by Eqs. (5)-(10) of Ref. Cowern et al., 1990 as

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Here we introduced the parameter $\nu$ to facilitate the proof that the profile Eq. (27) of Ref. Cowern et al., 1990 coincides with the 1D profile obtained from Eq. (19) by the inverse Fourier transform with $C_{0}(k)=1$ corresponding to the delta-function initial profile:

[TABLE]

Here the last exponential on the first line has been expanded in the Tailor series so by comparison with Eq. (27) we conclude that the inverse Fourier transforms on the second line should be equal to $\phi_{n>0}(\xi,1)$. To show this we first introduce the integrals

[TABLE]

and note that if $\phi_{1}(\xi,\nu)$ is known explicitly, other integrals can be computed recursively as

[TABLE]

Thus, we only need to show that $\phi_{n}(\xi,\nu)$ in Eq. (31) satisfy the recursion. To this end we first note that with the exponential factor being common to all terms in all functions, the equality in Eq. (34) will hold if it will be valid for every power of $|\xi|^{k}$ under the summation sign. Let us consider one such term in Eq. (31)

[TABLE]

When substituted in Eq. (34) it will contribute to $|\xi|^{k}$ term in $\phi_{n+1}^{(k)}$ through the derivative of its last factor with respect to $\nu$. The only other contribution from $\phi_{n}$ contributing into $|\xi|^{k}$ term in $\phi_{n+1}$ is $\phi_{n}^{(k-1)}$ differentiated with respect to $\nu$ in the exponential function. It is straightforward to check that these two contributions lead to the term $\phi_{n+1}^{(k)}$ as in Eq. (35) only with $n+1$ instead of $n$, as required.

Appendix B Diffusion near the boundary in 1D geometry

If only segregation near the surface is of interest and $\lambda\ll L$ (in the examples in the main text $\lambda=5-10$ nm and $L=200$ nm), at small $x$ the influence of the second surface of the film at $x=2L$ can be neglected and the system can be considered as infinite in the direction of $x>0$. In this case the diffusion kernel in the half-space $x,x_{0}\geq 0$ can be found by the method of images as

[TABLE]

It is straightforward to check that (36) satisfies both equation (12) and the boundary condition (13) at $x=0$.

Even in 1D systems with the simple stationary pair creation rate (23) the diffusion equation (11) cannot be solved analytically for all $t$. However, at the beginning of the evolution when the impurity distribution still has its simple constant initial form Eq. (24) the solution can be found to the first order in $t$ as the sum of the initial distribution and the integral on the r.h.s. of Eq. (11) multiplied by $t$.Cowern et al. (1990) Substituting the constant $C_{i0}$, Eqs. (23) and (36) in Eq. (11) one gets after simple integrations

[TABLE]

up to a small correction of $O(1/s)$.