Direct Kinetic Monte Carlo Simulations of Interdiffusion

Abstract

Kinetic Monte Carlo (KMC) simulations of the diffusion couple experiments were performed with the assumption that the vacancy composition in the system equilibrates much faster than the atomic configuration. Within this approach, the consistent atomistic simulation model with immediate vacancy equilibration mechanism was developed by incorporating a physical model of vacancy sources and sinks into the KMC algorithm. The Semi-Grand Canonical Monte Carlo (SGCMC) algorithm determined equilibrium vacancy composition and configuration in a system and, when implemented with the KMC code, generated on-line vacancy compositions locally equilibrated according to the atomic configuration in the sample. The values of the interdiffusion coefficients were determined by means of the Boltzmann-Matano formalism applied to the simulated composition profiles along the diffusion couple. The simulations also clearly reproduced the Kirkendall effect expected to appear in the simulated systems. Validity and reliability of the approach was assessed by comparing the results with the predictions of the Darken-Manning theory.

Introduction

Parameterization of diverse diffusion processes is based on the notion of a diffusion coefficient $D_i$ defined by the 1^st Fick law relating the flux ${\overrightarrow{j}}_i$ of a system component $i$ to the gradient $\overrightarrow{\mathrm{\nabla }}{C}_i$ of its composition:

1

$${\overrightarrow{j}}_i=-D_i·\overrightarrow{\mathrm{\nabla }}{C}_i$$

Equation 1 is written in the crystal frame of reference. When combined with the continuity equation $-\overrightarrow{\mathrm{\nabla }}·{\overrightarrow{j}}_i=\frac{\partial C_i}{\partial t}$ the 1st Fick law yields

2

$$\frac{\partial C_i}{\partial t}=\overrightarrow{\mathrm{\nabla }}·\left(D_i·\overrightarrow{\mathrm{\nabla }}{C}_i\right)$$

Equation 2 is referred to as the 2nd Fick law.

In general, the diffusion coefficient $D_i$ may vary with composition $C_i$, temperature, T, and pressure, P. Here, we consider only composition dependence $D_i=D_i\left(C_i\right)$.

The relation Eq (1) is used for description of two main diffusion quantities: (i) tracer diffusion coefficient $D_i^{\ast }$ for a process of a spontaneous mixing of atoms/molecules taking place at the absence of composition (or chemical potential) gradient; and (ii) intrinsic (chemical) diffusion coefficient $D_i^I$ for processes driven by gradients $\overrightarrow{\mathrm{\nabla }}{\mu }_i$ of chemical potentials.^[1] The well-known formula ${\mu }_i={\mu }_i^{\left(0\right)}-k_B·T·ln\left({\gamma }_i·C_i\right)$; where $k_B$ and ${\gamma }_i$ denote the Boltzmann constant and the activity coefficient of the component $i$, respectively, leads to the relationship between $D_i^I$ and $D_i^{\ast }$:

3

$$D_i^I=D_i^{\ast }·\left(1+\frac{\partial ln{\gamma }_i}{\partial ln{C}_i}\right),$$and to the definition of so-called thermodynamic factor ${\varphi }_i=1+\frac{\partial ln{\gamma }_i}{\partial ln{C}_i}$.

Among various atomistic mechanisms of diffusion^[2] of special importance is the “vacancy mechanism”, where atoms in a solid (e.g. crystalline) may move exclusively by jumping to neighbouring vacancies. The vacancy mechanism controls diffusion in most metallic systems (alloys, intermetallics) and was revealed due to the Kirkendall effect commonly observed in interdiffusion experiments.^[3]

It may be shown [1, p.171] that in a binary alloy with vacancies A-B-V where the vacancy concentration $C_V\ll C_{A\left(B\right)}$ the thermodynamic factors ${\varphi }_A$ and ${\varphi }_B$ are equal: ${\varphi }_A={\varphi }_B=\varphi$.

The fundamental relationship relating the tracer diffusion coefficient $D_{}^{\ast }$ defined by the 1^st Fick law (Eq. (1)) to the atomistic diffusion mechanism (e.g. the vacancy mechanism) is the Einstein-Smoluchowski formula:

4

$$D_{}^{\ast }=\underset{t\to 0}{lim}\frac{<r_i{\left(t\right)}^2>}{6t}$$where $<r_i{\left(t\right)}^2>$ denotes the mean square distance (MSD) travelled by an atom of type i within time $t$. (see e.g.^[4]).

The MSD may be evaluated by means of atomistic simulations or modelled within different atomistic approaches (see e.g.^[5]) involving frequencies $w_i$ of elementary atomic jumps (e.g. to vacancies), as well as the so-called correlation functions $f_i^I$ and $f_i$ , depending on the geometry of the jumps, and, also, on their mechanism.

Consequently, the tracer and intrinsic diffusion coefficients (Eq. (3)) can be written as

5

$$D_i^I=D_i^{\left(0\right)}·f_i^I·{\varphi }_i;D_i^{\ast }=D_i^{\left(0\right)}·f_i$$where $D_i^{\left(0\right)}$ represents a hypothetical coefficient of tracer diffusion running with purely random (uncorrelated) walk of a species $i$.

The tracer correlation function $f_i$ can be derived from the Einstein-Smoluchowski formula, Eq. 4, as follows:

6

$$f_i=\underset{t\to \infty }{lim}\frac{<r_i{\left(t\right)}^2>}{<r_i^{\left(\mathit{random}\right)}{\left(t\right)}^2>}$$where $<r_i^{\left(\mathit{random}\right)}{\left(t\right)}^2>$ denotes MSD corresponding to purely random (uncorrelated) walk of an average atom of species $i$.

Interdiffusion means a diffusional exchange of atoms/molecules across two different alloys that are in contact (a diffusion couple). The process is driven by the chemical potential gradient across the boundaries and if initially the two materials constituting the diffusion couple differ in compositions of the components, it usually acts towards homogenization of the diffusion couple (Fig. 1). Further description applies to the process running in a binary diffusion couple initially composed of two monatomic blocks A and B.

[See PDF for image]

Fig. 1

Scheme of the process of interdiffusion: diffusion couple A-B at $t=t_0$ (a) and at $t>t_0$ (b); grey, black and white points denote A-atoms, B-atoms and vacancies, respectively. A-atom composition profile $C_A\left(x\right)$ at $t=t_0$ and at $t>t_0$ is shown in (c)

The interface between the A- and B-components of the diffusion couple is represented by so-called Matano plane, where the net flux of A- and B-atoms equals zero^[1]. The Matano plane is located at the position $x_M$ (Fig. 2) defined by:

7

$$\underset{C_A=1}{\overset{0}{\int }}\left(x-x_M\right)d{c}_A=0$$

[See PDF for image]

Fig. 2

Position $x_M$ of the Matano plane. Eq. (6) implies that the areas $S_1$ and $S_2$ are equal

In the case of different diffusion rates of A- and B-atoms, the Matano plane moves with respect to the laboratory frame and in such a case is often called a Kirkendall plane. This effect was observed for the first time by Kirkendall and Smigelskas, in 1941 ^[3], and it is called now the Kirkendall effect.

The process of interdiffusion is parameterized by the composition-dependent interdiffusion coefficient $\tilde{D}\left(C_A\right)$ which fulfils Eq. (2) (re-written in the laboratory reference frame) at each point across the diffusion couple and parameterizes the flattening rate of the composition profile. Of interest is the determination of the coefficient $\tilde{D}$ which can be achieved in two ways: (i) directly, by analyzing the time dependent composition profile $C_A\left(x\right)$ across the diffusion couple (Boltzmann-Matano (B-M) technique)^[6,7]:

8

$$\tilde{D}\left(C_A^{\left(0\right)}\right)=-\frac{1}{2t}·\frac{{\int }_{C_A=1}^{C_A^{\left(0\right)}}\left[x\left(C_A\right)-x_M\right]·d{C}_A}{{\left(\frac{d{C}_A}{\mathit{dx}}\right)}_{C_A^{\left(0\right)}}}$$or (ii) by exploring the theoretical relationship between $\tilde{D}$ and the tracer diffusion coefficients $D_A^{\ast }$ and $D_B^{\ast }$ of the system components determined by means of independent experiments.

The latter approach is based on the atomistic Darken-Manning model of vacancy-mediated interdiffusion yielding (see e.g.^[1]):

9

$$\tilde{D}\left(C_A,C_B\right)=C_A·D_B^I+C_B·D_A^I$$where $D_A^I$ and $D_B^I$ are intrinsic diffusion coefficients of the components A and B (Eq. (3)).

In conditions of the vacancy-mediated interdiffusion, the relationship between $D_A^I$ and $D_B^I$ and the corresponding tracer diffusion coefficients $D_A^{\ast }$ and $D_B^{\ast }$ is, however, additionally affected by so-called ‘vacancy wind’ which appears when $D_A^{\ast }\ne D_B^{\ast }$. In such a case Eq. (5) holds with

10

$$f_A^I=f_{\mathit{AA}}-f_{\mathit{AB}}·C_A·C_B^{-1};f_B^I=f_{\mathit{BB}}-f_{\mathit{BA}}·C_B·C_A^{-1}$$where the ‘collective correlation factors’ $f_{\mathit{ij}}$ are equal to

11

$$f_{\mathit{ij}}=\underset{t\to \infty }{lim}\frac{<R_i\left(t\right)·R_j\left(t\right)>}{<R_i^{\left(\mathit{random}\right)}{\left(t\right)}^2>}$$where $R_i\left(t\right)$ denotes a sum of the displacements$r_i\left(t\right)$of all i-atoms performed within the time $t$ and averaging is performed over independent systems. Consequently, the collective correlation factors $f_{\mathit{ii}}$ (Eq.(11)) are not identical with the tracer correlation factors $f_i$ defined by Eq. (6), where averaging is performed over individual atoms.

The Darken-Manning formula (Eq. (9)) is often rewritten as

12

$$\tilde{D}\left(C_A,C_B\right)=\left(C_A·D_B^{\ast }+C_B·D_A^{\ast }\right)·\varphi ·S$$where $S$ denotes the ‘vacancy wind factor’:

13

$$S=\frac{C_A·D_B^{\ast }·s_B+C_B·D_A^{\ast }·s_A}{C_A·D_B^{\ast }+C_B·D_A^{\ast }}$$where $s_A=\frac{f_{\mathit{AA}}-C_A/C_B·f_{\mathit{AB}}}{f_A}$ and $s_B=\frac{f_{\mathit{BB}}-C_B/C_A·f_{\mathit{BA}}}{f_B}$.

Typically, $S$ is close to unity, however, $s_A$ and/or $s_B$ could be significantly different from unity.

Specific atomistic formalisms allowing the calculation of the above correlation factors have been developed:

1. the “assembly method” based on simulations,^{[8, 9–10]}

2. an analytical approach (referred to as MAA) yielding more accurate results but being applicable only to random alloys.^[11]

The aim of the present work was to elaborate a package of Monte Carlo (MC) algorithms for direct simulation of vacancy-mediated interdiffusion accompanied by the Kirkendall effect.

While the standard Kinetic Monte Carlo (KMC) algorithm is a natural tool for the simulation of vacancy-mediated atomic migration generating the time evolution of the composition profile $C_A\left(x\right)$ (Fig.1) to be analysed within the B-M formalism (Eq.(8)), a separate problem to be solved is vacancy composition and distribution $C_V\left(x\right)$ in the simulated diffusion couple.

As discussed in the previous papers (see e.g.^[12]), experimental and theoretical investigations reveal that the rate with which vacancy composition adjusts to the current chemical composition and configuration in a system – a process involving vacancy migration, vacancy creation and vacancy annihilation – is much higher than the rate of structural rearrangements controlled by atomic migration. Thus, it seemed reasonable to assume that the simulated interdiffusion would run with vacancies showing composition profile $C_V\left(x\right)$ adjusted – by local equilibration – to the chemical composition profile $C_A\left(x\right)$. Noteworthy, such conditions perfectly meet the basic assumption of the Darken-Manning model applied in the present study.

While the equilibrium chemical-composition-dependent vacancy composition $C_V\left(C_A\right)$ is determinable by means of Semi-Grand Canonical Monte Carlo (SGCMC) simulations,^[13] due to the generally non-linear $C_V\left(C_A\right)$ dependence, local equilibration of $C_V$ during isothermal interdiffusion requires that the number of vacancies present in the couple is constantly updated. This means, in turn, that the MC simulation of vacancy-mediated interdiffusion must be combined with a particular model for distribution of vacancy sources and sinks.

In view of all the above requirements, the software elaborated within the present work consists of three algorithms:

• SGCMC algorithm determining equilibrium vacancy composition $C_V\left(C_A\right)$;

• standard KMC algorithm simulating vacancy-mediated atomic migration;

• special MC algorithm combined with a model of vacancy sources and sinks generating on-line the profile $C_V\left(x\right)$ of vacancy composition equilibrated with respect to the current $C_A\left(x\right)$ profile.

The quality and correctness of the direct simulation of interdiffusion were evaluated by comparing the $\tilde{D}\left(C_A\right)$ dependence resulting from the B-M analysis of the simulated composition profiles $C_A\left(x\right)$ (Fig. 1) with $\tilde{D}\left(C_A\right)$ determined by means of the “assembly method” – i.e. considerations based on the Darken-Manning model (Eqs. (7) – (11)). This task gave rise to additional extended MC simulations of the homogeneous A-B binary systems aiming at the determination of the chemical composition dependences of the following equilibrium and kinetic parameters:

• the thermodynamic factors ${\varphi }_A\left(C_A\right)$ and ${\varphi }_B\left(C_A\right)$;

• the average frequencies $w_A\left(C_A\right)$ and $w_B\left(C_A\right)$ of A- and B-atom jumps to the vacancies residing at the nearest neighbour (nn) position;

• the tracer correlation factors $f_i\left(C_A\right),i=A,B$

• the tracer diffusion coefficients $D_A^{\ast }\left(C_A\right)$ and $D_B^{\ast }\left(C_A\right)$

• the collective correlation factors $f_{\mathit{ij}}\left(C_A\right),i,j=A,B$

Methodology

MC Simulation Techniques

Determination of Equilibrium Vacancy Composition and the Thermodynamic Factors as Functions of Chemical Composition

Equilibrium vacancy compositions in the systems were determined by means of Semi-Grand Canonical Monte Carlo (SGCMC) simulations.^[13] There, a bcc A-B-V lattice gas decomposing into a vacancy-rich and vacancy-poor phases is simulated. The vacancy-poor phase being in equilibrium with the vacancy-rich one is identified for the A-B-V system of interest and the composition of vacancies it contains is interpreted as the equilibrium vacancy composition. As a ternary system, the A-B-V lattice gas may decompose at a given temperature into phases showing a range of chemical compositions and thus, the method yields the equilibrium vacancy composition as a function of $C_A$. The above phase equilibrium automatically yields the (relative) values of the chemical potentials ${\mu }_A\left(C_A\right)-{\mu }_V\left(C_A\right)$ and ${\mu }_B\left(C_A\right)-{\mu }_V\left(C_A\right)$ and in view of ${\mu }_V\approx 0$ enables the calculation of the thermodynamic factors ${\varphi }_A\left(C_A\right)$ and ${\varphi }_B\left(C_A\right)$ defined by Eq. (2).

Simulation of Vacancy-Mediated Atomic Migration

The vacancy-mediated diffusion kinetics in the systems were simulated with the following approximation of the frequency of the jumps of a $p-$ atom residing on $i-$ lattice site to a vacancy residing on the nearest-neighbour (nn) $j-$ lattice site:

14

$${\mathrm{\Pi }}_{p,i\to j}={\mathrm{\Pi }}_0·exp\left[-\frac{E_{p,i\to j}^{+}-E_{p,ij}}{k_B·T}\right]$$where $E_{p,i\to j}^{+}$ denotes the saddle-point energy related to such jumps and $E_{p,ij}$ is the energy of the system with a p-type atom occupying the i-th lattice site and being a nearest neighbour (nn) of a vacancy residing on the j-th lattice site (see Fig. 3).

[See PDF for image]

Fig. 3

Variation of the system configuration energy due to a jump performed by a p-type atom initially occupying the i-th lattice site and moving to a nearest neighbour (nn) vacancy residing on the j-th lattice site

Similarly to the previous works (see e.g.^[14]), the saddle-point energy $E_{p,i\to j}^{+}$ was approximated by

15

$$E_{p,i\to j}^{+}=\frac{E_{p,ij}+E_{p,ji}}{2}+E_{\mathit{bar}}^{+}$$where $E_{p,ji}$ denotes the system energy after the $p$-type atom $i\to j$ jump – i.e. with a p-type atom occupying the j-th lattice site and being a nearest neighbour (nn) of a vacancy residing on the i-th lattice site. The migration barrier parameters $E_{\mathit{bar}}^{+}$ depended exclusively on the kinds of the jumping atoms.

The jump frequencies for A- and B-atoms (Eqs. 14 and 15) were implemented in a standard “residence-time” KMC algorithm which yielded: (i) the composition profile $C_A\left(x\right)$ across the simulated diffusion couple after a given number of MC-steps; (ii) the average frequencies $w_A$ and $w_B$ of the A- and B-atom jumps to nn vacancies and (iii) the MC-time dependence of mean-square-distances (MSD) $<r_A^2\left(t\right)>$ and $<r_B^2\left(t\right)>$ travelled by A- and B-atoms during the process of self-diffusion – the parameters used in the “assembly method” and MAA analysis.

Model of Vacancy Sources and Sinks

Basic Assumptions

The Kirkendall’s experiments have shown that neither surfaces, nor the grain boundaries serve as sources/sinks for vacancies, as the measured movement of the sample is independent of the specimen size and grain size ^[15]. This suggests the existence of many potential sources and sinks within a given grain in the form of dislocations. Bardeen and Herring^[16] showed that a mixed dislocation could operate as a source/sink by giving off or taking on a plane of vacancies if the system is supersaturated with vacancies up to 1%.

The model applied in the present paper is based on the mechanism of dislocation climbing. Two major phenomena are related to this process in real materials: (i) an atom-vacancy exchange and (ii) a local displacement due to the generated internal stress. According to the discussion presented in the previous section, use of the Darken-Manning model is now justified provided the rate of mechanical relaxation of stress due to dislocation sliding is assumed much higher than the rate of diffusion.

The effect of dislocation climbing and sliding can be illustrated by considering a 2D boundary case of a supersaturated system with vacancies, see Fig. 4. After many acts of the dislocation climbing up (removing vacancies) and following stress relaxation, the dislocation leaves the sample. As a result, an atomic plane is removed from the surface, and the sample shrinks by a size proportional to the lattice constant (Burgers vector). Practically, the infinitesimal internal displacements due to the annihilation of a single vacancy cause a displacement of the whole sample by $\frac{\mathit{dx}}{N_{\mathit{plane}}}$ , where $N_{\mathit{plane}}$ is a number of atoms on the single atomic plane and $\mathit{dx}$ is the change of the sample size when an entire plane of atoms is removed. In the opposite situation where vacancies are generated, the sample expands.

[See PDF for image]

Fig. 4

An example of the sample shrinking due to the annihilation of vacancies. Dashed line marks the initial size of the sample; $\mathit{dx}$ is a total change of the size when the dislocation is removed from the sample

An Algorithm for the Operation of Vacancy Sources

The main idea of the elaborated algorithm is to mimic the real process based on the dislocation climbing and sliding by a specific vacancy migration on a rigid lattice (with no dislocation). Fig. 5 shows a scheme of a single vacancy annihilation. The vacancy exchanges with the neighbouring atom and generates the first stage of a dislocation movement (climbing). In the next step, hypothetical local distortions are transported to the outer part of the sample (surface, inclusions, vacancy pore) by mimicking the dislocation sliding on the perfect lattice – Fig. 5b. This movement is realised by a series of consecutive vacancy- atom exchanges with the atoms located in the path of the possible dislocation movement. In this way, the sample shrinks by the volume decrement proportional to the volume of the single vacancy–Fig. 5(c). It is of course assumed that the volume of the vacancy is equal to the volume of the atom.

[See PDF for image]

Fig. 5

Scheme of a vacancy annihilation. (a) Change in the distortion in the real material due to the dislocation climbing up; (b) A real sample is mapped onto the perfect lattice. The path of a hypothetical dislocation sliding is shown; (c) The modelled sample after annihilation of the vacancy

Simulation of the dislocation climbing requires that: (i) the path of the dislocation movement is defined and (ii) a series of atom-vacancy exchanges in the direction of the path is executed. Definition of the path is based on the Astar algorithm^[17], which allows determination of the shortest path between the surface and a site in the bulk.

The context of interdiffusion requires in addition that the atomic/vacancy fluxes generated by the simulated process mimicking dislocation climbing and sliding do not affect the position of the Matano plane (see section 1). Fulfilment of this requirement is ascertained by an additional algorithm.

Detailed description of the algorithms elaborated and used in the present study is available in ref. [18].

KMC with On-Line Vacancy Equilibration Algorithm

Implementation of the model of vacancy sources and sinks required that the standard sample of a diffusion couple was completed with reservoirs of vacancies from which vacancies might be collected (and replaced by atoms), or where they might be deposited. In the present work, blocks of empty bcc unit cells were added to both sides of the couple (Fig. 6)

[See PDF for image]

Fig. 6

A simulated diffusion couple built of A-atoms (dark grey), B-atoms (grey) and vacancies (bright grey). The diffusion couple is wrapped between two blocks of empty bcc cells

The KMC simulations were stopped every $\mathrm{\Delta }n$ steps and unit cells of the diffusion couple were scanned layer by layer for current vacancy composition $C_V$ and chemical composition $C_A$. Depending on the sign of $\mathrm{\Delta }{C}_V=C_V-C_V^{\left(\mathit{eq}\right)}\left(C_A\right)$ the number $n=\mathrm{\Delta }{C}_V·N_{\mathit{layer}}$ (N_layer denoting the number of lattice sites in a single layer of unit cells) of either vacancies $\left(\mathrm{\Delta }{C}_V<0\right)$, or atoms $\left(\mathrm{\Delta }{C}_V>0\right)$ were chosen at random in the layer and were exchanged with atoms or vacancies, respectively, by means of the algorithm mimicking the dislocation movement (section 2.2.2.).

Simulated Systems

Hamiltonians and Atomic Migration Energies

Two binary bcc systems with vacancies, A_cB_(1-c)-V, were simulated: a system showing no long- range ordering (denoted as ‘DIS’) and a system with B2-ordering tendency (denoted as ‘ORD’). Similarly, as in our previous papers (see e.g.^[14]), vacancies were treated as a third component.

As only the ‘DIS’ system was treatable by means of the MAA theory, it served as a model system enabling definite assessment of the reliability of the used MC algorithms. Both the ‘DIS’ and ‘ORD’ systems were modelled with Ising Hamiltonians involving nearest-neighbor (nn) pair interaction energies V_pq whose values are displayed in Table 1.

Table 1. The nn pair interaction energies V_pq parameterizing the simulated A_cB_(1-c)-V systems

$V_{\mathit{pq}}$	$V_{\mathit{AB}}$	$V_{\mathit{AA}}$	$V_{\mathit{BB}}$	$V_{\mathit{AV}}$	$V_{\mathit{BV}}$	$V_{\mathit{VV}}$
‘DIS’ (eV)	−  0.11	−  0.11	− 0.11	0.0	0.0	0.0
‘ORD’ (eV)	− 0.10	− 0.05	− 0.05	+ 0.025	+ 0.025	0.0

Two sets of the applied migration barriers $E_{\mathit{bar}}^{+}\left(A\right)$ and $E_{\mathit{bar}}^{+}\left(B\right)$ (see Fig. 3) are shown in Table 2.

Table 2. Applied values of the migration barriers $E_{\mathit{bar}}^{+}\left(A\right)$ and $E_{\mathit{bar}}^{+}\left(B\right)$

System type/barrier type	‘SYM’		‘ASY’
	$E_{\mathit{bar}}^{+}\left(A\right)$	$E_{\mathit{bar}}^{+}\left(B\right)$	$E_{\mathit{bar}}^{+}\left(A\right)$	$E_{\mathit{bar}}^{+}\left(B\right)$
‘DIS’ (eV)	0.5	0.5	0.204	0.402
‘ORD’ (eV)	0.55	0.55	0.204	0.402

While the ’DIS’ system with the ’SYM’ barriers was ideal for testing the algorithm, the ’DIS’ system with the ’ASY’ barriers was applied to model the diffusion of a heavier isotope A* (here denoted as B) in the pure A solid. The ‘ORD’ systems with ’SYM’ and ‘ASY’ barriers represented a simple intermetallic and an intermetallic showing the Kirkedal effect, respectively.

As both the “assembly” and MAA formalisms involve the value of the ratio of the average frequencies $w_A$ and $w_B$ of the A- and B-atom jumps to nn vacancies, in order to make the planned analysis as clear as possible, the KMC simulations were performed at temperatures $T$ yielding the same value of $\frac{w_A}{w_B}=10$ for all the systems. This was achieved at $T=1000\text{K}$ for ‘DIS’ and at $T=1400\text{K}$ for ‘ORD’ at stoichiometric AB composition. In the case of the stoichiometric AB ‘ORD’ system the temperature $T=1400\text{K}$ corresponded to $\frac{T}{T_{O-D}}=0.8$ ( T_O-D denoting the “order-disorder” transition temperature) which meant that the simulated A_cB_(1-c)-V system exhibited B2-long-range order (B2 LRO) within a certain range of $c$.

Simulated Samples

Samples Generated for the SGCMC and KMC Simulations of Homogeneous A-B-V Systems

Following the procedure described in detail in ref. [19], the SGCMC simulations always started with stoichiometric AB supercells built of 30x30x30 bcc unit cells (i.e. containing $N_S=54000$ lattice sites) with 3D periodic boundary conditions. In order to avoid the problem of metastable long-lived antiphase domains (in the case of the ‘ORD’ system), the samples were initially built with maximum degree of B2 LRO – i.e. all of the A and B resided exclusively on the α- and β-sublattices, respectively (Fig.7).

[See PDF for image]

Fig. 7

Schematics of the B2-type superstructure: α-sublattice (empty circles); β-sublattice (solid circles)

The SGCMC simulations generated homogeneous samples of A_cB_(1-c)-V with equilibrium vacancy composition in which vacancy-mediated self-diffusion was simulated by means of the KMC algorithm.

Samples Generated for the Simulations of Interdiffusion (Diffusion Couples)

Diffusion couples were composed of two initially monatomic A- and B-blocks built of 45 × 25 × 25 bcc unit cells ($N_S=56000$ lattice sites) each and of two blocks of 10 layers of empty (vacant) unit cells attached to both ends of the couple (Fig. 6). The periodic boundary conditions

were imposed only in the directions perpendicular to the composition gradient.

Results

Parameters Yielded by MC Simulations of Homogeneous A_cB_(1-c)-V Systems

Equilibrium Parameters

Equilibration of the simulated systems was made by saturation of the monitored MC time dependence of their internal energy. Equilibrium values of the considered parameters were determined by averaging over 5000 observations separated by 200000 Monte Carlo steps from the saturated region.

The effect of B2 LRO in the ‘ORD’ A_cB_(1-c) clearly shows up in the c-dependence of the equilibrium compositions of A- and B-anti-site atoms (Fig. 8).

[See PDF for image]

Fig. 8

Chemical composition dependences of the compositions $C_A^{\left(\beta \right)}$ of A-anti-sites (solid circles) and $C_B^{\left(\alpha \right)}$ of B-anti-sites (open circles) (a) in the ‘DIS’ systems at $T=1000\text{K}$; and (b) in the ‘ORD’ systems at $T=1400\text{K}$ (b). $C_A^{\left(\beta \right)}=\frac{2·N_A^{\left(\beta \right)}}{N_S}$, $C_B^{\left(\alpha \right)}=\frac{2·N_B^{\left(\alpha \right)}}{N_S}$ ($N_A^{\left(\beta \right)}$ and $N_B^{\left(\alpha \right)}$ are the numbers of A-atoms residing on the $\beta$-sublattice and B-atoms residing on the $\alpha$-sublattice, respectively)

Loss of the linearity of $C_A^{\left(\beta \right)}\left(c\right)$ and $C_B^{\left(\alpha \right)}\left(c\right)$ in the ‘ORD’ systems for $0.3<c<0.7$ reflects the B2 LRO occurring in this range of $c$ at $T=1400\text{K}$.

Figure 9 shows the composition dependence of the equilibrium vacancy composition $C_V^{\left(\mathit{eq}\right)}\left(c\right)$ in both simulated systems. While lack of the composition dependence of $C_V^{\left(\mathit{eq}\right)}$ in the ‘DIS’ system with $V_{\mathit{AB}}=V_{\mathit{AA}}=V_{\mathit{BB}}$ verifies the correctness of the SGCMC simulations, strong composition dependence caused by $V_{\mathit{AB}}\ne V_{\mathit{AA}}$ is observed in ‘ORD’ system: $C_V^{\left(\mathit{eq}\right)}$ reaches minimum value at the stoichiometric composition $c=\frac{1}{2}$, which is in agreement with our previous results (see e.g.^[19]).

[See PDF for image]

Fig. 9

Chemical composition dependence of the equilibrium vacancy composition $C_V^{\left(e,q\right)}$: (a) in the ‘DIS’ systems at $T=1000\text{K}$; and (b) in the ‘ORD’ systems at $T=1400\text{K}$

The final result of the SGCMS simulations – the composition dependence of the thermodynamic factors ${\varphi }_A$ and ${\varphi }_B$, is shown in Fig. 10. In the case of the ‘DIS’ system, which is an ideal solution both ${\varphi }_A$ and ${\varphi }_B$ are composition independent and approximately equal to unity. This again verifies the correctness of the approach. In the ‘ORD’ system ${\varphi }_A\left(c\right)$ and ${\varphi }_B\left(c\right)$ show distinct maxima within $0.3<c<0.7$, where B2 LRO means different occupations of the neighbouring lattice sites and thus, an increase of the gradients $\overrightarrow{\mathrm{\nabla }}{\mu }_{A\left(B\right)}$ of the chemical potentials of diffusing atoms. The difference between the evaluated ${\varphi }_A$ and ${\varphi }_B$ is negligible and this agrees with theoretical predictions (see section I).

[See PDF for image]

Fig. 10

Chemical composition dependence of the thermodynamic factors ${\varphi }_A$ (open circles) and ${\varphi }_B$ (open squares) in: (a) the ‘DIS’ systems $T=1000\text{K}$; and (b) the ‘ORD’ systems at $T=1400\text{K}$

Kinetic Parameters Obtained by KMC Simulations

Figure 11, 12 and 13 show the composition dependences of the average A- and B-atom jump frequencies $w_A$ and $w_B$, as well as their ratios $w_A/w_B$ in the ‘DIS-ASY’, ‘ORD-SYM’ and ‘ORD-ASY’ A_cB_(1-c) systems with vacancies. Obviously, no composition dependence of the jump frequencies $w_A$ and $w_B$ is observed in both ‘DIS’ systems. However, while $w_A\approx w_B$ held in the ‘DIS-SYM’, in the ’DIS-ASY’ the ratio $w_A/w_B$ was close to 10 (Fig. 11), as earlier assumed.

[See PDF for image]

Fig. 11

Chemical composition dependence of the ratio $w_A/w_B$ of A- and B-atom average jump frequencies in the ‘DIS-ASY’ systems at $T=1000\text{K}$

[See PDF for image]

Fig. 12

Chemical composition dependences of: (a) average A-atom and B-atom jump frequencies $w_A$ (solid circles) and $w_B$ (empty circles), respectively; and (b) the ratio $w_A/w_B$, in the ‘ORD-SYM’ systems at $T=1400\text{K}$

[See PDF for image]

Fig. 13

Chemical composition dependences of: (a) average A-atom jump frequency $w_A$; (b) average B-atom jump frequency $w_B$; and (c) the ratio $w_A/w_B$ (c) in the ‘ORD-ASY’ systems at $T=1400\text{K}$

Strong composition dependences of $w_A$, $w_B$ and $w_A/w_B$ were observed in ‘ORD-SYM’ and ‘ORD-ASY’ systems (Fig. 12 and 13). In both cases the jump frequencies showed minimum values at the stoichiometric composition $c=1/2$ which clearly corroborated with minimum vacancy composition (Fig. 9b). The desired value of $\frac{w_A}{w_B}=10$ was achieved in ‘DIS-ASY’ only for $c=1/2$.

The composition dependences of the tracer diffusion coefficients $D_A^{\ast }\left(C_A\right)$ and $D_B^{\ast }\left(C_A\right)$ resulting from KMC simulations of ‘DIS-ASY’, ‘ORD-SYM’ and ‘ORD-ASY’ systems with equilibrium vacancy compositions are shown in Fig. 14, 15 and 16.

[See PDF for image]

Fig. 14

Chemical composition dependences of tracer diffusion coefficients: (a) $D_A^{\ast }$; and (b) $D_B^{\ast }$, in the ‘DIS-ASY’ A_cB_(1-c) systems at $T=1000\text{K}$

[See PDF for image]

Fig. 15

(a) Chemical composition dependences of tracer diffusion coefficients $D_A^{\ast }$ (solid circles) and $D_B^{\ast }$ (empty circles) in the ‘ORD-SYM’ A_cB_(1-c) systems at $T=1400\text{K}$. (b) Chemical composition dependence of the ratio $D_A^{\ast }/D_B^{\ast }$

[See PDF for image]

Fig. 16

Chemical composition dependences of tracer diffusion coefficients: (a) $D_A^{\ast }$; and (b) $D_B^{\ast }$, in the ‘ORD-ASY’ A_cB_(1-c) systems at $T=1400\text{K}$. (c) Chemical composition dependence of the ratio $D_A^{\ast }/D_B^{\ast }$

The tracer diffusion coefficients $D_A^{\ast }\left(C_A\right)$ and $D_B^{\ast }\left(C_A\right)$ were evaluated as

16

$$D_{A\left(B\right)}^{\ast }=\underset{t\to \infty }{lim}\left[\frac{1}{6N_{A\left(B\right)}}·\frac{\partial }{\partial t}\left(\sum _{N_{A\left(B\right)}}{r}_{A\left(B\right)}^2\left(t\right)\right)\right]$$

Figure 14 shows the self-diffusion coefficients for the ’DIS-ASY’ system, where the well marked difference between the composition dependences of $D_A^{\ast }$ and $D_B^{\ast }$ (despite almost constant $w_A/w_B$ (Fig. 11)) reflects correlation effects (see Fig. 17 below).

[See PDF for image]

Fig. 17

Chemical composition dependences of the tracer correlation factors: (a) $f_A$; and (b) $f_B$, in the ‘DIS-ASY’ system at $T=1000\text{K}$

The composition dependences of $D_A^{\ast }$ and $D_B^{\ast }$ in the ‘ORD-SYM’ and ‘ORD-ASY’ systems (Fig. 15 and 16) reflect the influence of the B2 LRO which generally suppresses self-diffusion (see the minima of $D_A^{\ast }$ and $D_B^{\ast }$ at $c=1/2$) – the effect is discussed in ref. [19]. Asymmetry of the migration barriers $E_{\mathit{bar}}^{+}\left(A\right)$ and $E_{\mathit{bar}}^{+}\left(B\right)$ generally increased the ratio $D_A^{\ast }/D_B^{\ast }$ without, however, affecting qualitative features of $D_A^{\ast }\left(c\right)$ and $D_B^{\ast }\left(c\right)$.

The composition dependences of the tracer correlation factors (Eq. (6)) following from the KMC simulations of the ‘DIS’ systems at $T=1000\text{K}$ and of the ‘ORD’ systems at $T=1400\text{K}$ are displayed in Figs. 17,18 and 19. Low correlation factor indicates that high fraction of atomic jumps to vacancies were reversed. In the case of ‘DIS-ASY” the relationship $f_A<f_B$ (Fig.17) apparently follows from $E_{\mathit{bar}}^{+}\left(A\right)<E_{\mathit{bar}}^{+}\left(B\right)$ causing more A-atom jumps being statistically reversed than the B-atom ones.

[See PDF for image]

Fig. 18

Chemical composition dependences of the tracer correlation factor, $f_A$ (solid circles) and $f_B$ (empty circles) in the ‘ORD-SYM’ system at $T=1400\text{K}$

[See PDF for image]

Fig. 19

Chemical composition dependences of the tracer correlation factors: (a) $f_A$; and (b) $f_B;$ at $T=1400\text{K}$ in the ‘ORD-ASY’ system

In the case of the ’ORD’ systems the composition dependences of $f_A$ and $f_B$ show well-marked anomalies at the limits of the range of LRO at $T=1400\text{K}$ (Figs. 18 and 19). The lower values of $f_A$ and $f_B$ in the structure with B2 LRO – at the minimum values at $c=1/2$ marking the highest degree of LRO – resulted from the definite enhancement of reversing the A- and B-atom jumps to anti-site positions (see e.g.^[19] for more detailed discussion).

Figs. 20 and 21 show the composition dependences of the collective correlation factors $f_{\mathit{AA}}$, $f_{\mathit{BB}}$, $f_{\mathit{AB}}$ and $f_{\mathit{BA}}$ defined by Eq. (11). In all cases, it was observed that $\underset{C_i\to 0}{lim}{f}_{\mathit{ii}}=f_i$, $\underset{C_i\to 1}{lim}{f}_{\mathit{ii}}=1$ and $\underset{C_i\to 1}{lim}{f}_{\mathit{ij}}=0$. This is consistent with the definition in Eq. (11) which for a single particle case reduces to the Eq. (6) and with the fact that a minority component only weakly affects the diffusion of the majority atoms.

[See PDF for image]

Fig. 20

Chemical composition dependences of the collective correlation factors $f_{\mathit{AA}}$ (solid circles), $f_{\mathit{BB}}$ (empty circles), $f_{\mathit{AB}}$ (solid squares) and $f_{\mathit{BA}}$ (empty squares) in the: (a) ‘DIS-SYM’ and (b) ‘DIS-ASY’ systems at $T=1000\text{K}$

[See PDF for image]

Fig. 21

Chemical composition dependences of the collective correlation factors $f_{\mathit{AA}}$ (solid circles), $f_{\mathit{BB}}$ (empty circles), $f_{\mathit{AB}}$ (solid squares) and $f_{\mathit{BA}}$ (empty squares) in the: (a) ‘ORD-SYM’ and (b) ‘ORD-ASY’ systems at $T=1400\text{K}$

In general, it was observed that $f_{\mathit{ii}}$ increased and $f_{\mathit{ij}}$ decreased with increasing $C_i$ composition. In the case of the ordered systems the composition dependences of the collective correlation factors showed clear minima at $c=0.5$ (Fig.21)) which was again an effect of the predominance of reversal jumps in the diffusion process in B2 superstructure.

Interdiffusion Coefficients Calculated Within the “Assembly” Formalism

In the case of the ’DIS-SYM’ system the interdiffusion coefficient was constant and equal to

the self-diffusion coefficients $\tilde{D}=D_A^{\ast }=D_B^{\ast }$ – as expected from the Eq. (9) applied to this case.

Figures 22a, 23a and 24a show the composition dependencies of the interdiffusion coefficients $\tilde{D}$ calculated for ’DIS-ASY’, ’ORD-SYM’ and ’ORD-ASY’ systems according to Eqs. (4), (9) and (10).

[See PDF for image]

Fig. 22

Chemical composition dependences of: (a) the interdiffusion coefficient $\tilde{D}$; and (b) the vacancy wind factor $S$; calculated for the ‘DIS-ASY’ system within the “assembly” formalism. Dotted line marks the $\tilde{D}\left(c\right)$ curve yielded by the MAA theory^[11]

[See PDF for image]

Fig. 23

Chemical composition dependences of: (a) the interdiffusion coefficient $\tilde{D}$; and (b) the vacancy wind factor $S$; both are calculated for the ‘ORD-SYM’ system within the “assembly” formalism

[See PDF for image]

Fig. 24

Chemical composition dependences of: (a) the interdiffusion coefficient $\tilde{D}$; and (b) the vacancy wind factor $S$; both are calculated for the ‘ORD-ASY’ system within the “assembly” formalism

In the case of the ’DIS-ASY’ system the vacancy wind factor S > 1 (Fig.22b) indicates a strong influence of the correlation effects resulting from $D_A^{\ast }>D_B^{\ast }$ holding in the whole composition range. Thanks to the random character of this system, it was possible to compare the “assembly” results with those yielded by the Moleko-Allnatt-Allnatt theory (MAA)^[11] (see the dotted line in Fig. 22a). Perfect agreement between both $\tilde{D}\left(c\right)$ dependences verifies correct implementation of the software developed in the present study.

In the case of ordering systems, the diffusion of the atomic components and, consequently, the interdiffusion in the ordered region ($0.3<c<0.7$, Figs. 23a and 24a) was slowed down with respect to the diffusion in a disordered structure.

Despite a strong maximum of the thermodynamic factor $\varphi ={\varphi }_A={\varphi }_B$ at $c=1/2$ (Fig. 10b)

the evaluated interdiffusion coefficient $\tilde{D}\left(c\right)$ shows minimum indicating its prevailing dependence on the vacancy wind factor S (Figs 23b and 24b). The strong minimum of $S\left(c\right)$ at $c=1/2$ is, in turn, strictly correlated with a change of sign of $D_A^{\ast }-D_B^{\ast }$ (Fig.15a) which generates additional flux of vacancies on both sides of the initial contact plane.

In the case of the ’ORD-ASY’ system, the correlation effects influence the interdiffusion process in a more complex way by slowing it down in the regions reached by the slower component (c < 0.5) and by accelerating it in the region reached by the faster component (c > 0.5). This effect was correlated, in turn, with a sudden increase of the ratio $D_A^{\ast }/D_B^{\ast }$ at c > 0.5 (Fig. 16c).

Results of Direct Simulations of Interdiffusion

General Assumptions

Each ‘simulation experiment’ covered 100 independent simulations of the generated diffusion couple (section 3.2.2.) by means of the algorithm described in section 2.2. The composition profile, fluxes of the constituents and the general configuration of the simulated couple were monitored in the simulations.

The profiles resulting from separate simulations were merged and sorted according to the MC

time. Subsequently, averaging over time and space was performed by means of the ‘weighted

moving average procedure’.^[20]

The B-M analysis of the composition profiles $C_A\left(x\right)$ corresponding to the MC-times $t_i$ (Fig. 1, Eq. (6)) yielded the chemical composition dependence $\tilde{D}\left(C_A\right)$ of the interdiffusion coefficient related to $t_i$.

The final result of the analysis is presented in Fig.25.

[See PDF for image]

Fig. 25

MC-time evolution of the interdiffusion coefficient $\tilde{D}\left(C_A,=,0.5\right)$ in the ’DIS-SYM’ system at $T=1000\text{K}$

Despite initially evolving with MC-time the evaluated interdiffusion coefficient $\tilde{D}$ always saturated and thus, the final value was obtained by averaging over the level of saturation resulting from the given simulation. The developed software is available on-line: https://github.com/ufsowa/tools_direct.git

Simulation Without Vacancy Sources and Sinks

Only the ’DIS-ASY’ diffusion couple was simulated without vacancy sources and sinks. The couple composed of two equal blocks of A- and B-atoms, respectively, contained equilibrium number of vacancies determined by means of SGCMC simulations of the stoichiometric AB system (see section 2.1.1). The vacancies were initially uniformly distributed along the couple and their total number was constant all over the simulation.

The composition dependence of the simulated interdiffusion coefficient is presented in Fig. 26.

[See PDF for image]

Fig. 26

Chemical composition dependence of the interdiffusion coefficient $\tilde{D}$ in the ’DIS-ASY’ yielded by B-M analysis of the composition profiles simulated with constant vacancy composition at $T=1000\text{K}$: raw results of these simulations (open circles) and averaged values of $\tilde{D}$ (solid circles)

The difference with respect to the corresponding theoretical result presented in the Fig. 22 is significant. The value of $\tilde{D}$ (e.g. the rate of mixing) in the simulated diffusion couple was the highest in the region with $C_A=1$ and was decreasing monotonically to the lowest value in the region with $C_A=0$. This effect was a consequence of the distribution of vacancies across the diffusion zone of the couple shown in Fig. 27.

[See PDF for image]

Fig. 27

Distributions of A-atoms (solid line) and vacancies (cloud of points) along the simulated ‘DIS-ASY’ diffusion couple at $T=1000\text{K}$. Equilibrium vacancy composition (independent of the chemical composition (Fig.9a)) is marked by a dotted line

Due to the asymmetry of the A- and B-atom jump frequencies $w_A=10·w_B$, a net flux of vacancies occurred and drained the vacancies from the B-rich side to the A-rich side of the couple generating highly non-equilibrium vacancy distribution in the diffusion zone (Fig.27). Low vacancy composition on the B-rich side obviously damped the vacancy-mediated interdiffusion. The considerable discrepancy between the “assembly” and direct simulation results obviously follows from the violation of the basic Darken-Manning assumption of the equilibrium vacancy distribution in the system.

Simulation with the Vacancy Equilibration Algorithm

The equilibration of vacancy composition (see section 2.2) was performed every 1000 KMC steps which ensured that the correlation effects were saturated^[21] and equilibrium composition and distribution of vacancies was approximately achieved (Fig. 28a). The created and annihilated atoms and vacancies were monitored and analyzed. Figure 28b shows that vacancies were brought in and taken out mostly within the diffusion zone (i.e. near the initial contact plane) and that the algorithm worked correctly–i.e. any locally occurring shortage or excess of vacancies was properly compensated.

[See PDF for image]

Fig. 28

(a) Vacancy distribution developed in the diffusion couple of ’DIS-ASY’ system simulated with the equilibration mechanism at $T=1000\text{K}$; (b) vacancies generated and annihilated in the sample after MCtime = 4500 s. Empty circles represent the deficit/excess of vacancies $\mathrm{\Delta }{N}_V=N_V-N_V^{\left(e,q\right)}$ generated by interdiffusion; solid circles show the number $N_V$ of vacancies brought in ($N_V>0$) or taken out ($N_V<0$) by the equilibration algorithm

Chemical composition dependencies $\tilde{D}\left(c\right)$ of the interdiffusion coefficient in ‘DIS-ASY’ system, obtained by means of the ‘assembly’ method and the B-M analysis of the simulated $C_A\left(x\right)$ profiles are shown in Fig. 29.

[See PDF for image]

Fig. 29

Interdiffusion coefficient for ’DIS-ASY’ system simulated with vacancy equilibration mechanism at $T=1000\text{K}$ (solid circles). Open circles show results obtained with ‘assembly’ method (see Fig. 22a). Similarly, as in the Fig. 22a dotted line marks the $\tilde{D}\left(c\right)$ curve yielded by the MAA theory^[11]

Implementation of the vacancy equilibration algorithm resulted in a definite change of the B-M results (with respect to those displayed in Fig.26) which now showed a reasonable agreement with the ‘assembly’ ones. The still observed small discrepancy might arise from some deviation of the composition and configuration of vacancies from equilibrium, which can be seen in Fig. 28a revealing a small deviation from equilibrium close to the initial contact plane (x = 0).

Direct simulation of vacancy mediated interdiffusion in the ’ORD-ASY’ system yielded much smaller composition gradient which generated large inaccuracy of the B-M analysis. Nevertheless, the resulting c dependence $\tilde{D}\left(c\right)$ of the interdiffusion coefficient displayed in Fig. 30 shows a reasonable agreement with the corresponding ‘assembly’ result (Fig. 24a).

[See PDF for image]

Fig. 30

Chemical composition dependence of the interdiffusion coefficient $\tilde{D}$ in ’ORD-ASY’ system simulated with vacancy equilibration mechanism at T = 1400 "K" (solid circles) and yielded by the ‘assembly’ calculations (open squares)

The Kirkendall Shift

From the description of the equilibration algorithm it can be deduced that the implemented dislocation-mediated mechanism for vacancy sinks/sources can lead to a shrinking/expansion of the system when vacancies are annihilated or created, respectively. In the case when vacancies are created on one side and removed from the system on the opposite side, a movement of the sample should be observed. Indeed, this effect was clearly visible in the case of ’DIS-ASY’ system. The position $x_K$ of the Kirkendall (Matano) plane was found by monitoring a plane showing constant atomic composition $C_A=0.5$. Its displacement $\mathrm{\Delta }{x}_K$ with respect to the laboratory frame is displayed vs. MC-time in Fig. 31. The linearity of $\mathrm{\Delta }{x}_K\left[{\left(\text{MC - time}\right)}^2\right]$ is consistent with experimental observations^[3] and the predictions of Darken approximation.

[See PDF for image]

Fig. 31

MC-time dependence of the displacement $\mathrm{\Delta }{x}_K$ of the Kirkendall plane simulated in the ’DIS-ASY’ system at $T=1000\text{K}$

Summary

It was demonstrated that interdiffusion coefficient $\tilde{D}$ parameterizing the process of vacancy-mediated atomic migration in a diffusion couple may be evaluated by atomistic simulations. The principal goal was to elaborate a method for direct simulation of a diffusion couple and the determination of $\tilde{D}$ by means of the B-M analysis of the composition profile along the diffusion path. The procedure should be considered reliable and corresponding to a real process thanks to its implementation with a model of vacancy sinks and sources providing continuous local equilibration of vacancy composition and configuration. A simple model of a dislocation acting as a sink/source of vacancies within the rigid lattice approximation was proposed. Reliability of this new ’direct’ technique was tested by comparing the results with those yielded by the traditional analytical approaches: the MMA analytical method and the standard methodology based on the Darken-Manning theory and KMC simulations of self-diffusion in a homogeneous system. The simulated bcc systems were modelled with the Ising-type Hamiltonian.

The compatibility of the results obtained by means of both methods showed that: (i) the new model of vacancy sources/sinks was developed correctly and (ii) the Darken-Manning theory of diffusion and the related methods of simulations are valid in the limit of ${\mu }_V\approx 0$.

It was shown that the direct method could be a good solution in cases where thermodynamic quantities are studied in a dynamically evolving system.

It is planned to study the thermodynamic quantities off equilibrium by considering component activities, as well as the augmented distribution of vacancy sources/sinks. It will also be possible to study diffusion at the conditions of complex gradients and system geometries. This is important for newly elaborated theories of diffusion–e.g. in case of multi-component systems or diffusion with vacancy sources/sinks.^{[22, 23–24]}

Acknowledgments

Financial assistance from Polish National Science Center [Grant No. 2015/16/T/ST3/00501], Polish Ministry of Science and High Education [Grant No. 3135/7. PR/2014/2], European Community’s Seventh Framework Programme (FP7 PEOPLE-2013-IRSES) [Grant No. EC-GA 612552] and the Endeavour Fellowship [Grant No. ERF-RDDH-5049-2016] and the Australian Research Council Discovery Grants [Grant No. ARC DP200101969] schemes are greatly acknowledged.

Data Availability

The raw and processed data required to reproduce these findings are available to download from https://github.com/ufsowa/tools_direct.git

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

1. Mehrer, H. Diffusion in Solids: Fundamentals, Methods, Materials, Diffusion–Controlled Processes; 2007; New York, Springer: 168. [DOI: https://dx.doi.org/10.1007/978-3-540-71488-0] 1138.76050

2. Manning, JR. Diffusion Kinetics for Atoms in Crystals; 1968; Princeton, N.J., Van Nostrand: 8.0426.58016

3. Kirkendall, E; Smigelskas, A. Zinc Diffusion in Alpha Brass. Trans. Metall. Soc. AIME; 1947; 171, pp. 130-142.

4. Islam, MA. Einstein-Smoluchowski Diffusion Equation: A Discussion. Phys. Scr.; 2004; 70, 2–3 pp. 120-125.2004PhyS..70.120I [DOI: https://dx.doi.org/10.1088/0031-8949/70/2-3/008] 1139.82325

5. Bakker, H. Murch, GE; Nowick, AS. Tracer Diffusion in Concentrated Alloys. Diffusion in Crystalline Solids; 1984; Orlando, Academic Press: pp. 189-256. [DOI: https://dx.doi.org/10.1016/B978-0-12-522662-2.50009-1] 0706.73025

6. Boltzmann, L. Zur Integration Der Diffusionsgleichung Bei Variabeln Diffusions Coefficienten. Ann. Phys.; 1894; 289, 13 pp. 959-964. [DOI: https://dx.doi.org/10.1002/andp.18942891315] 25.1553.02 in German

7. Matano, C. On the Relation Between the Diffusion-Coefficients and Concentrations of Solid Metals (the nickel-copper system). Jpn. J. Phys.; 1933; 8, 3 pp. 109-113.0445.35063

8. Qin, Z; Murch, GE. Computer Simulation of Chemical Diffusion in a Binary Alloy with an Equilibrium Concentration of Vacancies. Philos. Mag. A; 1995; 71, 2 pp. 323-332.1995PMagA.71.323Q [DOI: https://dx.doi.org/10.1080/01418619508244359] 0834.11050

9. I.V. Belova and G. E. Murch, A Theory of Diffusion in the Ordered Alloy III. Chemical Diffusion, Philos. Mag. A, 1996, 73(6), p 1699–1713.

10. Van der Ven, A; Yu, H; Ceder, G; Thornton, K. Vacancy Mediated Substitutional Diffusion in Binary Crystalline Solids. Prog. Mater. Sci.; 2010; 55, 2 pp. 61-105. [DOI: https://dx.doi.org/10.1016/j.pmatsci.2009.08.001] 1056.74057

11. Moleko, LK; Allnatt, AR; Allnatt, EL. A Self–Consistent Theory of Matter Transport in a Random Lattice Gas and Some Simulation Results. Philos. Mag. A; 1989; 59, 1 pp. 141-160.1989PMagA.59.141M [DOI: https://dx.doi.org/10.1080/01418618908220335] 1006.82510

12. R. Kozubski, D. Kmieć, E. Partyka, and M. Danielewski, ‘Order-order’ Kinetics in Ni50.5Al49.5 Single Crystal, Intermetallics, 2003, 11, p 897–905.

13. Biborski, A; Zosiak, L; Kozubski, R; Sot, R; Pierron-Bohnes, V. Semi-Grand Canonical Monte Carlo Simulation of Ternary Bcc Lattice-Gas Decomposition: Vacancy Formation Correlated with B2 Atomic Ordering in A-B Intermetallics. Intermetallics; 2010; 18, pp. 2343-2352. [DOI: https://dx.doi.org/10.1016/j.intermet.2010.08.007]

14. Sowa, P; Kozubski, R; Biborski, A; Levchenko, EV; Evteev, AV; Belova, IV; Murch, GE; Pierron-Bohnes, V. Self Diffusion and ‘Order-Order’ Kinetics in B2-Ordering AB Binary Systems with a Tendency for Triple Defect Formation: Monte Carlo Simulation. Philos. Mag.; 2013; 93, pp. 1987-1998.2013PMag..93.1987S [DOI: https://dx.doi.org/10.1080/14786435.2012.742591]

15. P. Shewmon, Diffusion in Solids. The Minerals, Metals & Materials Series, Springer, 2016, p 132–145.

16. J. Bardeen and C. Herring, Diffusion in Alloys and the Kirkendall Effect, Atom Movements, J. H. Hollomon, Ed., Cleveland, OH, American Society for Metals, 1951, p. 87–111.

17. Philibert, J. Diffusion and stresses. Defect Diffus. Forum; 1996; 129, pp. 3-8. [DOI: https://dx.doi.org/10.4028/www.scientific.net/DDF.129-130.3] 1066.14005

18. P. Sowa. PhD Thesis, Jagiellonian University 2017, Krakow, Poland. https:/ruj.uj.edu.pl/xmlui/handle/item/310841

19. Betlej, J; Sowa, P; Kozubski, R; Murch, GE; Belova, IV. Self-Diffusion in a Triple-Defect A-B Binary System: Monte Carlo simulation. Comput. Mater. Sci.; 2020; 172, [DOI: https://dx.doi.org/10.1016/j.commatsci.2019.109316] 109316.

20. Kouya, S; Miwa, S. A Method of Scr Improvement with a Two-Dimensional Weighted Moving Average Filter. Electron. Commun. Japan; 2003; 86, pp. 58-67. [DOI: https://dx.doi.org/10.1002/ecja.1172] 1323.65082

21. Belova, IV; Murch, GE. Collective Diffusion in the Binary Random Alloy. Philos. Mag. A; 2000; 80, pp. 599-607.2000PMagA.80.599B [DOI: https://dx.doi.org/10.1080/01418610008212070] 1054.53019

22. S. Kornienko, A. Gusak, and G. Lutsenko, Nonequilibrium vacancies in nanosystems, in Diffusion and Stresses, vol. 264, Trans Tech Publications, 2007, p 109–116

23. Svoboda, J; Fischer, F; Fratzl, P. Diffusion and Creep in Multi-Component Alloys with Non-Ideal Sources and Sinks for Vacancies. Acta Mat.; 2006; 54, pp. 3043-3053.2006AcMat.54.3043S [DOI: https://dx.doi.org/10.1016/j.actamat.2006.02.041] 1426.74136

24. Allnatt, AR; Paul, TR; Belova, IV; Murch, GE. A High Accuracy Diffusion Kinetics Formalism for Random Multicomponent Alloys: Application to High Entropy Alloys. Phil. Mag.; 2016; 96, pp. 2969-2985.2016PMag..96.2969A [DOI: https://dx.doi.org/10.1080/14786435.2016.1219785]