Self-diffusion in a triple-defect A-B binary system : Monte Carlo simulation

(1)

Contents lists available atScienceDirect

Computational Materials Science

journal homepage:www.elsevier.com/locate/commatsci

Self-di ﬀusion in a triple-defect A-B binary system: Monte Carlo simulation

J. Betlej

^a

, P. Sowa

^a

, R. Kozubski

^a,⁎

, G.E. Murch

^b

, I.V. Belova

^b

aM. Smoluchowski Institute of Physics, Jagiellonian University, Lojasiewicza 11, 30-348 Krakow, Poland

bCentre for Mass and Thermal Transport in Engineering Materials, School of Engineering, The University of Newcastle, Callaghan 2308, Australia

A R T I C L E I N F O

Keywords:

A. intermetallics

B. diﬀusion, thermodynamic properties D. point defects

E. Monte Carlo simulation, defects: theory

A B S T R A C T

In this comprehensive and detailed study, vacancy-mediated self-diffusion of A- and B-elements in triple-defect B2-ordered ASB1-Sbinaries is simulated by means of a kinetic Monte Carlo (KMC) algorithm involving atomic jumps to nearest-neighbour (nn) and next-nearest-neighbour (nnn) vacancies. The systems are modelled with an Ising-type Hamiltonian with nn and nnn pair interactions completed with migration barriers dependent on local configurations. Self-diffusion is simulated at equilibrium and temperature-dependent vacancy concentrations are generated by means of a Semi Grand Canonical MC (SGCMC) code. The KMC simulations reproduced the phenomena observed experimentally in Ni-Al intermetallics being typical representatives of the 'triple-defect' binaries. In particular, they yielded the characteristic‘V’-shapes of the isothermal concentration dependencies of A- and B-atom diffusivities, as well as the strong enhancement of the B-atom diffusivity in B-rich systems. The atomistic origins of the phenomenon, as well as other features of the simulated self-diffusion such as temperature and composition dependences of tracer correlation factors and activation energies are analyzed in depth in terms of a number of nanoscopic parameters that are able to be tuned and monitored exclusively with atomistic simulations. The roles of equilibrium and kinetic factors in the generation of the observed features are clearly distinguished and elucidated.

1. Introduction

The notion of the‘triple defect’ was introduced by Wasilewski[1]

who originally deﬁned it as a complex of a single A- or B-antisite defect and two nn vacancies in a stoichiometric A-50 at% B system with the B2 superstructure (Fig. 1). Generation of ‘triple defects’ stems from a substantial diﬀerence between the formation energies for A- and B- antisite defects whose consequence is that the system disorders (e.g.

due to increasing temperature) by preferentially creating the antisites with lower formation energy. Such a phenomenon is called‘triple-defect disordering’ (TDD). It should be noted that in general TDD determines only statistics of the generated defects which may occur without the generation of compact‘triple defects’ as deﬁned by Wasi- lewski[1].

The tendency for TDD implies: (i) a large diﬀerence between the A- and B-antisite concentrations; (ii) a large diﬀerence between the concentrations of vacancies residing on α- and β-sublattices (the ‘home’

sublattices of A and B atoms); (iii) an increase of vacancy concentration with decreasing degree of chemical long-range order– i.e. with an increasing concentration of antisite defects. In the extreme case of the exclusive generation of A-antisites, their concentration is equal to one half of the vacancy concentration – i.e. the vacancy concentration

strongly increases with decreasing degree of chemical order. In non- stoichiometric binaries the tendency for TDD – i.e. lower formation energy for A-antisites, means that while A-antisites compensate for the deficit of B atoms in A-rich systems, the B-atoms in B-rich systems re- main on the β-sublattice and the departure from stoichiometry is compensated by‘structural’ α-vacancies. The process of TDD should be contrasted from the so called triple-defect mechanism of diffusion which means atomic migration via specifically correlated atomic jumps mediated by vacancy-pairs[2].

The topic of self or tracer diﬀusion in stoichiometric and non-stoichiometric A-B B2 intermetallics with the tendency for TDD has been widely investigated. Except for the basic interest in the physical aspects of the phenomenon, the studies were taken up due to technological attraction of the most common TDD alloys (NiAl, FeAl, CoAl,…). As aluminides, they are successfully applied in industrial manufacturing and in coatings exposed to high temperature in aggressive and corrosive environments.

While the number of theoretical and computational studies is fairly large, experimental works are relatively rare. The main reason for this lies in the major diﬃculties posed by experimental tracer-diﬀusion methods that require radioactive isotopes of the constituents. Many TDD alloys have aluminum as the second component which

https://doi.org/10.1016/j.commatsci.2019.109316

Received 14 May 2019; Received in revised form 19 September 2019; Accepted 26 September 2019

⁎Corresponding author.

E-mail address:rafal.kozubski@uj.edu.pl(R. Kozubski).

Available online 14 October 2019

T

(2)

unfortunately lacks suitable radioactive isotopes. The direct tracer-diffusion experiments concerned, therefore the transition-metal components– mostly Ni, whereas the tracer diffusion coefficient of Al has been mainly deduced (approximately) from interdiffusion experiments, the transition-metal self-diffusion coefficient and the thermodynamic factor.

In 1949, Smoluchowski and Burgess[3]measured the tracer diffusion coefficient of Co in NiAl. Co is known for its general tendency to substitute for Ni on the Ni sublattice. Radioactive Co was plated on the NiAl sample and the decrease of activity due to the penetration of co- balt into material was recorded at 1150 °C. The overall shape of the graph showing the self-diffusion coefficient of Co vs. concentration of Ni perfectly resembles one obtained in a more recent experiment by Frank et al.[4].

In 1971, Hancock and MacDonnell[5]measured the tracer diffusion coefficient of the radioactive isotope⁶³Ni in NiSAl1-Spolycrystals over the temperature range of 1270–1600 K. The lowest value of diffusion coefficient was observed for the stoichiometric alloy (S = 0.5), with a relatively high activation energy of 3.2 eV. It appeared that even a small excess of Al led to a rapid decrease of the activation energy to 1.84 eV for S = 0.483. Almost symmetrically, for Ni-rich material, tracer diffusion was also faster than in the stoichiometric system, with a gradual increase of diffusion coefficient and a lowering of the activation energy towards 2.25 eV for S=0.58.

These results were critically evaluated by Frank et al.[4]. The most significant improvement with respect to the aforementioned study was the almost exclusive use of NiSAl1-S monocrystals. Accordingly, the impact of grain boundaries was minimized leading to diffusion coefficients that were smaller by about an order of magnitude. The crucial result of this experiment was that while the Ni-tracer diffusivity sys- tematically rose withSno significant change of Ni tracer diffusivity was observed in the range of0.468<S<0.5.

Much less information is available concerning Al diffusion. In 1975 Lutze-Birk and Jacobi[6]measured¹¹⁴In tracer diffusion in NiAl. Being in the same group as Al in the periodic table,¹¹⁴In replaces Al on the Al sublattice. As a function of chemical composition the¹¹⁴In-diffusivity showed the characteristic‘V’-shape with the minimum value around the stoichiometric composition.

Indirect estimation of the tracer diffusion coefficients in Ni-Al from interdiffusion experiments has recently been performed by Paul et al.

[7]and Minamino et al.[8].Fig. 2shows the results obtained by Paul et al. which suggests that when traced in a logarithmic scale versus concentration both Ni and Al tracer diﬀusivities show again the ‘V’- shape with minima around the stoichiometric composition Ni-50at.

%Al. Another important feature is the intersection of theDNiandDAl

isotherms atS<0.5.

Predicting rapid growth of the tracer diﬀusivities of Ni and Al with an increase of Ni and Al content respectively, the results of Paul and Minamino are in good qualitative agreement with those of Frank et al.

[4]and Hancock and MacDonnell[5]. In addition, the absolute values

of the tracer diﬀusion coeﬃcients calculated for Ni and In[6]are very close to the quantities directly measured in the vicinity ofS=0.5.

However, the symmetrical growth of the Ni diﬀusivity with decreasing S (the‘V’-shape) is in clear contrast with the most reliable results of Frank et al.[4].

The above experimental results have been widely analyzed in terms of the activation energy of diﬀusion and possible mechanisms of atomic migration responsible for this energy. Such mechanisms operating in Ni-Al intermetallics at temperatures at which experimental studies are performed are determined by a very high degree of the B2 long-range order (LRO) maintained in these systems up to the melting point of about 1900 K[9]. Krachler et al.[10]remarked on the impact of short- range chemical order in Ni-Al resulting in the curved shape of the Ar- rhenius plots of the measured tracer diﬀusivities [4]. Mishin et al.

[11–13] performed extensive studies of the energetics of point defect complexes and migration barriers for various diffusion mechanisms in Ni-Al modelled with interatomic potentials determined within the embedded atom method (EAM). Soule De Bas and Farkas[14]further extended that research by considering complex sequences of 10 and 14 atomic jumps. More recently, Chen et al.[15]and Yu et al.[16]analyzed the atomic migration barriers in Ni-Al applying either angle-dependent interactions or new EAM potentialsfitted to the experimental data. Marino and Carter proposed a more direct computational approach based solely on density functional theory (DFT) and in a series of works[17,18]evaluated not only the migration barriers and activation energies, but also the diffusion coefficients related to particular mechanisms proposed for Ni tracer diffusion in NiAl.

In 2011 Evteev et al.[19]published an interesting paper showing the results on Ni- and Al-self diffusion in NiAl simulated directly by means of Molecular Dynamics. The process was simulated in a layer limited by [1 1 0]-oriented free surfaces through which vacancies en- tered the system from outside and reached an equilibrium concentration. The simulations were performed at a temperature close to the melting point and yielded a Ni-diffusivity ca. 2.5 times higher than the Al diffusivity.

An almost complete computational study of diffusion in NiAl was performed by Xu and Van der Ven[20–22]who combined ab initio energy calculations with configurational thermodynamics by means of the Cluster Expansion method. The developed model of Ni-Al covered the equilibrium vacancy thermodynamics with, however, a priori assumptions concerning their preferential residence on particular sublattices in the B2 superstructure. Equilibrium thermodynamics of the system including vacancy concentrations was determined by means of the Semi Grand Canonical Monte Carlo (SGCMC) method assuming a Fig. 1. Scheme of B2-type superstructure: α-sublattice (unfilled circles); β-

sublattice (ﬁlled circles).

Fig.2. Reduced tracer diffusion coefficients of Ni and Al in NiSAl_1-S intermetallics at 1000 °C deduced from interdiffusion experiments. After Paul et al.

[7]

(3)

zero value of the chemical potential of vacancies. Separately, migration barriers were calculated in the DFT formalism for different jump types and local configurations. Consequently, it was possible to use kinetic Monte Carlo (KMC) to simulate diffusion processes in a system with a well-defined defect concentration and the degree of chemical order.

Thefinal results concerning the isothermal concentration dependence of Ni- and Al-diffusivities were in qualitative agreement with the experimental study of Frank et al.[4]reproducing the growth of Ni-diffusivity with increasing Ni-concentration in Ni-rich binaries. In agreement with the experimental works of Paul and Minamino [7,8]the same behaviour of Al-diffusivity has been observed. Much less attention has been paid to the Al-rich systems. The results shown in the work concern only one composition (C_Ni≈0.47) and may suggest that the diffusivity growth for both Ni and Al in the Al-rich region is much weaker than that reported experimentally[5–7]. Exploring the tracer- diffusion computations, Xu and Van der Ven evaluated the interdiffusion coefficient for Ni-Al which, however, decreased with growing Al content belowS=0.5. Although this clear contradiction with experiment could be attributed to the polycrystalline character of samples analyzed in the works [7,8], the authors suggested that it rather resulted from incorrect assumptions concerning the equilibrium vacancy concentration.

The present work aims at the determination and detailed analysis of the impact of the tendency for TDD– deﬁned at the beginning of this section, on self-diﬀusion of the components in B2-ordering A-B binaries.

Such studies have been taken up in the past (see e.g.[23]) focusing on the effect of interatomic interactions on the features of diffusion of system component atoms. By adapting in the model relationships between the atomic-jump migration energies yielded by ab-initio calculations dedicated to Ni-Al[22], the present study refers specifically to this system. The choice facilitates also the assessment of the simulation findings as the related experimental results with which they might be compared concern almost exclusively Ni-Al. The presented simulations address, therefore, vacancy-mediated atomic migration processes in a B2 superstructure of a TDD system loosely resembling Ni-Al.

By applying a straight forward Ising-type model it is possible to clearly demonstrate the strict correlation between the equilibrium thermodynamics of the system (equilibrium configurations of atoms and vacancies) and the kinetics of self-diffusion. Systems were simulated that represented uniformly a wide range of compositions both in the A-rich and B-rich side of the AB stoichiometry. The approach pro- vides a deep understanding of the diffusion phenomenon which is crucial for any effective development of material technologies.

The paper is organized as follows: The methodology of the study is described in Section 2 clearly pointing at the equilibrium and non- equilibrium (kinetic) aspects of the modelled phenomenon. The model of the simulated A-B binary showing the tendency for TDD and resembling the Ni-Al intermetallic system is presented inSection 3. The results of the study shown in detail inSection 4are then widely discussed inSection 5. The main conclusions are listed inSection 6.

2. Methodology

2.1. General remarks

The methodology of the reported study covers two aspects:

•

The determination of the temperature and composition dependence of the equilibrium atomic and point-defect conﬁgurations in the system.

•

The determination of the temperature and composition dependence of self-diﬀusivities and tracer correlation factors of the system components, as well as their activation energies.

In both cases, Monte Carlo (MC) simulations were performed.

Supercells were composed of 25 × 25 × 25 unit cells of the B2

superstructure (Fig. 1) – i.e. containing N=31250 lattice sites be- longing to equi-numerousα- and β-sublattices and populated with NA

A-atoms, NBB-atoms and NVvacancies. 3D periodic boundary condition (PBC) were imposed upon the supercells.

2.2. Model for equilibrium conﬁguration of the system

Of interest are the atomic conﬁgurations of a binary A-B system with vacancies. The conﬁgurations cover both the distribution of atoms over lattice sites and vacancy concentration and are parameterized by means of the following quantities:

•

Concentrations of atoms and vacancies on particularμ-sublattices:

= = =

C N

N (X A, B, V;μ α β, )

X

μ X

μ

( ) ( )

(1) where N_X^{( )}^μ denotes the number of X species residing onμ-sublattice

•

Total concentrations of atoms and vacancies:

= +

CX CXα C

Xβ

( ) ( )

(2)

•

Indicator of the chemical composition:

= +

+ + +

S N N

N N N N

Aα Aβ Aα

Aβ Bα

Bβ

( ) ( )

( ) ( ) ( ) ( )

(3)

•

Long-range order parameters for atoms and vacancies:

= −

+ = −

+

η C C

C C η C C

C C

, ,

A Aα

A β

Aα

Aβ B B

β Bα Bβ

Bα V Vα

Vβ Vα

Vβ

( ) ( )

(4)

•

Pair-correlations (short-range order parameters) for atoms and vacancies:

=

∙

C N

XY

μν XYμν

X totμν

( ) ( )

( )

(5) where N_XY⁽^μν⁾(X,Y = A,B,V;μ,ν = α,β) denotes the number ofX−Y nn or nnn pairs with X- and Y-species residing onμ−and −ν type sublattice sites, respectively;N_tot⁽^μν⁾= ∑_{X Y}_, N_XY⁽^μν⁾

It is important to note that because N is a sum of the numbers of atoms and vacancies, aﬁxed value of S does not imply constant values of the component concentrations CX deﬁned by Eq.(1).

In the equilibrium state corresponding to given external conditions particular configurations appear in the system with a specific probability distribution yielding average values of the above parameters interpreted as observables. By means of the Monte Carlo simulations it is possible tofind the equilibrium state of the system by generating a set of configurations showing the equilibrium probability distribution.

The present study was based on the Schapink model for the equilibrium configuration of a multicomponent system with vacancies[24], whose simple version was previously applied by one of the authors [25,26]. In this approach, a lattice gas A-B-V is treated as a regular ternary system – i.e. vacancies are treated strictly as an additional chemical component. The crucial property of the lattice gas (and also the condition for the applicability of the model) is that it shows a miscibility gap with a critical temperatureT_Cbelow which it decom- poses into two phases: one withCV ≪1and another (unrealistic) one withC_V≈1. Then, the basic assumption of the model is that the lattice- gas phase withC_V≪1being in equilibrium with the one withC_V ≈1is identified with the binary A-B crystal in equilibrium – i.e. the crystal with an equilibrium atomic configuration and equilibrium vacancy concentration.

(4)

2.3. Search for phase equilibria in the A-B-V lattice gas

Following the idea of Binder et al.[27]equilibrium compositions and configurations of the A-B-V lattice gas were determined at fixed temperatures T for arbitrary values of the chemical potentials μ_X (X = A,B,V). The procedure aimed atfinding their values μX

eq

( )yielding two solutions: one withC_V≪1and another withC_V≈1. Similar to our previous papers (see e.g. [28]) the lattice gas was examined using a standard algorithm of Semi Grand Canonical Monte Carlo (SGCMC) simulations where due to aﬁxed value of N the system is parameterized by two independent relative chemical potentials deﬁned in the present paper as:

= −

μ μ μ

Δ_AV _A _V (6)

= −

μ μ μ

Δ _BV _B _V

The SGCMC algorithm works in the following scheme:

(i) Random choice of a single lattice site occupied by a species ‘X’

(X = A,B,V);

(ii) Random choice of a species type‘Y’ (Y = A,B,V);

(iii) Replacement of the species ‘X’ by the species ‘Y’ with the Metropolis probability:

= ⎧

⎨⎩

⎡

⎣⎢− − − ⎤

⎦⎥

⎫

⎬⎭

→ →

exp E μ μ

Π min 1, Δ k T(Δ Δ )

X Y

X Y X Y

B (7)

whereΔE_X→_Y denotes the change of the system conﬁgurational energy due to theX→Y replacement. The quantity kBdenotes the Boltzmann constant.ΔE_X→_Y is evaluated within a particular model of the system implemented with the simulations and depends on the current composition and conﬁguration of the lattice gas.

(iv) Return to step (i).

Two series of SGCMC runs were performed at each temperature: in series 1 the simulations started with a perfect B2-ordered supercell with

= = = = =

C_A^{( )}^α C_B^{( )}^β ¹;C_A^β C_B^α C_V 0

2

( ) ( ) , whereas in the simulations of series 2 the supercell was initially empty (CA=CB=0;CV=1).

The SGCMC simulations run at temperatures belowTCyielded ty- pical C_V(Δμ_A, Δμ_B)isotherms as shown inFig. 3. The almost cliff-like discontinuity of the CV(ΔμA, ΔμB)surface reflected the coexistence of the vacancy-rich and vacancy-poor phases. The effect showed well- marked hysteresis and thus, exact evaluation ofΔμ_A^{( )}êq andΔμ_B^{( )}êq (the white line on theΔμ_A−Δμ_B plane) required an application of some further technique (see e.g.[28,29]). In the present work the technique of thermodynamic integration was chosen (see [28]for detailed de- scription and references).

2.4. Model of the vacancy-mediated atomic migration

Vacancy-mediated self-diffusion of A- and B-atoms was simulated at constant temperaturesTby means of the standard Residence-Time KMC algorithm[30]in samples withfixed chemical compositions (S) and equilibrium vacancy concentrationsCV corresponding toS andT and determined by the SGCMC runs. The initial atomic and vacancy configurations of the samples were generated by former SGCMC runs – i.e.

the initial values of CX( )μ

and CXV(μν)

(X = A,B,V;μ,ν = α,β) were close to the equilibrium (average) ones.

In view of the fact that in most of the previous papers devoted to the modelling of diﬀusion mechanisms in B2-ordering intermetallics, in particular, in Ni-Al, atomic jumps to both nn and nnn vacancies were considered, the same was implemented in the KMC algorithm applied in the present study.

The probability for an atom X (X = A,B) to jump from the initial i lattice site to a vacancy residing on nn or nnn j lattice site (Fig. 4) is given by:

= × ⎡

⎣

⎢− ⎤

⎦

⎥

→

exp E →

Π_{X i} _j Π k T^{X i} ^j

m

B

, 0

, ( )

(8) where:Π0is a pre-exponential factor whose value depends on the jump- attempt frequency of the X-atom and thus is, in general, a function of temperature and the type of jumping atom. The KMC-time increment of:

∑

= ×⎧

⎨⎩

⎡

⎣

⎢− ⎤

⎦

⎥

⎫

⎬⎭

→

−

t τ exp E

Δ k T

Xij

X im j B

:

( ) 1

(9) is assigned to each executed atomic jump.

E_{X i}^{( )}^m_,→_j is the migration barrier for the considered jump and according toFig. 4:

= −

→ +→

E_{X i}( )^m, _j E_{X i} _j EX ij

: : (10)

Within the model used in the present paper:

= +

→ +

+ → → +

E E E

E X

2 ( )

X i j

X i j X j i bar ,

( ) : :

(11) where the value of E_bar⁺ ( )X depends exclusively on the type X of jumping atom.

Consequently,

= −

→ +

→ → +

E E E

E X

2 ( )

X im j X j i X i j bar ,

( ) : :

(12) As was mentioned in our previous works (see e.g.[31,32]) such parameterization ofE_{X i}^{( )}^m_,_→_j partially accounts for its dependence on a local conﬁguration around the atom-vacancy pair.

Fig. 3. Typical CV(Δ , ΔμA μB)isotherm with a facet showing the coexistence of the vacancy-rich and vacancy-poor phases. The white line marks the positions ofΔμ_A^{( )}^eq andΔμ_B^{( )}^eq.

Fig. 4. Scheme of the energy parameterisation of a jump of a X-type atom (blue solid circle) residing on a lattice site i to a vacancy (black open square) residing on a lattice site j.

(5)

While justiﬁcation of the negligence of the temperature dependence ofΠ0(Eq. (1)was discussed earlier[32], almost equal values of the jump-attempt frequencies reported for Ni and Al-atoms (see e.g.

[22,33]) make it reasonable to assume in the present work a constant value ofΠ0equal to unity both for A- and B-atoms.

Finally, it should be noted that fulﬁlment of the detailed balance condition[34]by the KMC algorithm guaranteed conservation of the equilibrium conﬁgurations of the samples (i.e. maintenance of constant average values of CX( )μ

and CXV(μν)

) all over the KMC simulation runs performed atﬁxed temperatures.

2.5. Evaluation and analysis of diffusivities and correlation factors The self-diffusion coefficients DX for X-atoms (X = A, B) were evaluated from the standard Einstein-Smoluchowski relationship (see e.g.[35]):

= ⎡

⎣ 〈 〉⎤

⎦

→∞

D lim d

dt R t 1

6 ( )

X t X2

(13) where 〈R_X²( )t〉 denotes the monitored mean-square-distance (MSD) travelled by X-atoms (X = A, B) within the MC-time t– i.e. the value of R_X²( )t averaged over all X-atoms in the sample.

Analysis of the evaluated diﬀusivities in terms of the dynamics of atomic jumps to vacancies was done within the model of Bakker and colleagues[36]now extended upon atomic jumps to both nn and nnn vacancies. Expression of 〈R_X²( )t 〉in Eq. (13) in terms of elementary atomic jumps leads to

= ⎡

⎣⎢

〈 〉 × + 〈 〉 × ⎤

⎦⎥×

→∞

D lim n t a n t a

t f

{ ( ) ( ) }

X 6

t

Xnn

nn Xnnn

nnn X

corr

( ) 2 ( ) 2

( )

(14) where 〈n_X⁽ⁿⁿ⁾( )t〉and 〈n_X⁽ⁿⁿⁿ⁾( )t 〉denotes the average numbers of nn and nnn jumps performed by an X-atom within the MC-timet;a_nnand annn

denote the distances of the nn and nnn jumps, respectively; f_X⁽^corr⁾ denotes the tracer correlation factor given by[37]:

= ⎡

⎣⎢

〈 〉

〈 〉 × + 〈 〉 ×

⎤

⎦⎥

→∞

f lim R t

n t a n t a

( )

( ) ( )

X corr

t

X Xnn

nn Xnnn

nnn

( ) 2

( ) 2 ( ) 2

(15) The problem is conveniently parameterized with average atomic- jump frequenciesw_μ^{( )}^X_→_νdeﬁned as average numbers of jumps performed by one X-atom fromμ-sublattice sites to vacancies residing on ν-sublattice sites (μ,ν = α, β) within a unit KMC-time. Values ofw_μ^{( )}^X→_νare directly determined by counting the particular X-atomic jumps executed within aﬁxed number of KMC steps and by dividing the number of these jumps by the related KMC time interval and the numberN_X of X-atoms present in the supercell. Within the microscopic model[36]

they are expressed in terms of the atom-vacancy pair correlations CXV(μν)

(Eq. (5)) and the migration energies associated with the elementary atomic jumps (Eq.(8)):

= × × × ⎡

⎣

⎢−

〈 〉 ⎤

⎦

⎥

→

w S T z C S T exp E → S T

( , ) Π ( , ) k T( , )

μXν

μν XVμν X μm ν

B

( ) 0 ( ) ( ),

(16) where:

zμνdenotes the number ofν-sublattice sites being nn ( ≠μ ν) or nnn ( =μ ν) of aμ-sublattice site;

〈E_{X μ}^{( )}^m_,_→_ν( ,S T)〉 denotes the average over the migration barriers E_{X i}^{( )}^m_,→_jassociated with the jumps yieldingw_μ^{( )}^X→_ν( ,S T).

Due to the steady-state character of the simulated self-diﬀusion – i.e.

conservation of the average atomic conﬁguration (i.e. of the values of

〈C_X^{( )}^μ〉and〈C_XV⁽^μν⁾〉) guaranteed by the KMC algorithmw_μ^{( )}^X_→_ν _{fulﬁls the} detailed balance condition:

→ = →

w_μ^{( )}^X_ν w_ν^{( )}^X_μ (17)

Eqs.(16) and (17)yield, therefore, a link between the system energetics 〈{E_{X μ}^{( )}^m_, _→_ν〉} and the conﬁguration parameters {C_XV⁽^μν⁾}. Conse- quently, they also determine the steady-state atomic conﬁguration of the system at a temperatureT.

Eq.(17)implies that

〈 〉

= 〈 〉

= +

→ → →

n

t w n

t w w

2 ; [ ]

Xnn

αXβ Xnnn

αXα βXβ ( )

( )

( ) ( )

(18) Combination of Eqs.(14) and (18)yields:

= × × _→ × + _→ + _→ × ×

D 1 w a w w a f

6 {2 [ ] }

X α β

X

nn αXα β β

X

nnn X

corr

( ) 2 ( ) ( ) 2 ( )

(19) Eq.(19)makes it possible to demonstrate contributions of particular atomic jumps to the observed diffusion coefficients and thus to analyze in such terms the features of the effectively observed self-diffusion.

3. Model of the simulated system 3.1. Hamiltonian

Applied was an Ising-type model of the B2-ordering binary A-B system with vacancies (and, consequently, of the A-B-V ternary lattice gas) with nearest-neighbour (nn) V{ _XY⁽¹⁾} and next-nearest-neighbour (nnn) V{ _XY⁽²⁾}pair-interactions between atoms and vacancies: (X,Y = A, B,V). It should be noted that because of the varying composition of the system the SGCMC algorithm involves the total conﬁgurational energy (not only the energy of mixing) and therefore, separate evaluation of all the individual pair interaction parameters (not only of the‘ordering energies’W_XY^{( )}^j =2V_XY^{( )}^j −V_XX^{( )}^j −V_YY^{( )}^j) was required.

Evaluation of the V{ _XY⁽¹⁾} and V{ _XY⁽²⁾} parameters was based on the following criteria to be fulﬁlled by the modelled ternary A-B-V lattice gas:

(i) Ternary miscibility gap with a non-zero critical temperatureTC. (ii) B2-ordering of the vacancy-poor lattice-gas phase at temperatures

below the order-disorder temperatureT_{O D}− :T_{O D}₋ <T_C.

(iii) Tendency for TDD in the vacancy-poor phase – i.e. preferential formation of A-antisite defects. With reference to the earlier remarks (seeSection 1) the tendency for TDD was parameterized by the‘triple defect indicator’TDIdeﬁned in Ref.[36]asTDI= ^N

N Aβ

V ( )

. In the stoichiometric AB binary with the tendency for TDD

≈

TDI 1/2 should hold through aﬁnite temperature range.

As fulfilment of the above criteria determines only the relationships between V{ _XY⁽¹⁾}and V{ _XY⁽²⁾} assignment of particular values of the pair potentials required an arbitrary evaluation ofTO D− in the stoichiometric system withN_A=N_B. It should befirmly stressed that by no means did the latter affect meaningful results of the study, which in most cases are presented with relative (reduced) parameters.

The preliminary search for the proper values of V{ _XY⁽¹⁾}and V{ _XY⁽²⁾}was done by scanning their space and analytically checking the above criteria within the Bragg-Williams approximation (see[38]). As a starting point, the values of V{ _XY⁽¹⁾}found in Ref.[38]were used. A relationship

= − ×

V_XY⁽²⁾ 0.5 V_XY⁽¹⁾was chosen as an arbitrary assumption. Further ad- justment was performed by checking the equilibrium atomic conﬁg- urations generated by SGCMC simulations for the fulﬁlment of the criteria listed above. It must be emphasized that no calculations are known by the authors that accounted for interactions with vacancies and no strict reference to literature data was possible.

Theﬁnal values of V{ _XY⁽¹⁾}and V{ _XY⁽²⁾}used in the study are displayed in Table 1.

3.2. Migration barriers (saddle-point energies)

For the sake of the studies of atomic migration the extended Ising

(6)

model was completed with four parameters responsible for atomic migration: E_{nn bar}⁺_, ( )A,E_{nn bar}⁺_, ( )B, E_{nnn bar}⁺ _, ( )A andE_{nnn bar}⁺ _, ( )B (seeFig. 4 and Eq.(12)).E_{nn bar}⁺_, ( )X andE_{nnn bar}⁺ _, ( )X denote the parameters associated with X-atom jumps to nn and nnn vacancies, respectively.

While the relation E_{nn bar}⁺_, ( )A <E_{nn bar}⁺_, ( )B followed from the previous work [31], reference to the values ofE_{X i}^{( )}^m_,_→_j determined by ab initio calculations of Xu and Van der Ven[22](Table 2) suggested that E_{nnn bar}+ _, ( )A >E_{nnn bar}⁺ _, ( )B.

Although the evaluation of theE_{nn bar}⁺_, ( )X andE_{nnn bar}⁺ _, ( )X parameters was achieved by respecting the above criteria, the choice of the particular values was arbitrary with a lower limit yielded by the obvious condition of E_{X i}^{( )}^m_,_→_j>0. Table 3 displays the ﬁnal values of the E_{nn bar}+_, ( )X andE_{nnn bar}⁺ _, ( )X parameters used in the present KMC simulations.

4. Results

4.1. General remarks

Because of the arbitrary evaluation of the energetic parameters of the simulated A-B system, presentation of the simulation results in terms of absolute quantities was generally avoided. The values of particular parameters were, therefore, normalized to selected characteristic values as listed inTable 4.

4.2. Adequacy and eﬀectiveness of the Schapink model

Fig. 5shows an example of the isothermal S-dependence of the re- lative chemical potentials μΔ _AV^{( )}êqand μΔ _BV^{( )}êq determined both by directly analyzing the‘map’ of Fig. 3 [28]and by calculating their values in equilibrium configurations generated by standard canonical MC simulations[29].

Almost perfect agreement between the values obtained by both methods indicates the correctness of the present approach.

The SGCMC simulations of the A-B-V lattice gas revealed a miscibility gap, whose severalS=const. sections are shown inFig. 6.

The equilibrium B2-ordered phases with equilibrium vacancy concentrations were analyzed in the range of0.3<S<0.58 and showed S-dependent ‘order-disorder’ transition temperaturesTO D− never ex- ceeding the value ofT_C and reaching the maximum valueT_{O D}⁽^max₋ ⁾ for

≈

S 0.4(Fig. 7).

4.3. Point defect concentrations and TDD tendency of the system The SGCMC-generated equilibrium conﬁgurations of the system were analyzed in terms of the T- and S-dependence of point defect concentrations C_X^{( )}^μ (Eq.(1)), as well as of atom-vacancy pair correla- tions C_XV⁽^μν⁾(Eq.(5)). Arrhenius plots of the parameters were linear in a wide range of temperatures; however, showed well marked curvatures in the vicinity of the‘order-disorder’ transition points (seeFig. 8for the plots ofC_V, CV( )μ

and CXV(μν)). The eﬀect, obviously following from the increasing temperature dependence of the degree of LRO and SRO, was especially enhanced in the case ofα-vacancies in B-rich binaries.

As brieﬂy described inSection 1, the departure from stoichiometry of a B2 binary A-B system with a tendency for TDD is compensated by structuralα-vacancies in A-poor binaries ( <S 0.5) and by structural A- antisites in A-rich ones ( >S 0.5). The S-dependence of the concentrations of the structural point-defects atT→0K (i.e. at the absence of thermally activated defects) follows from the balance ∑_{X ν}_, N_X^{( )}^ν =N (X = A,B,V;ν = α,β) which yields:

(i) For structural vacancies ( <S 0.5;N_A^{( )}^β =N_B^{( )}^α =0):

Table 1

Values of nn (V_XY⁽¹⁾) and nnn (V_XY⁽²⁾) pair interaction parameters used in the study.

X-Y V_XY⁽¹⁾(eV) V_XY⁽²⁾(eV)

A-A –0.12 +0.06

B-B –0.05 +0.02

V-V 0 0

A-B –0.125 +0.062

A-V +0.038 –0.01

B-V –0.001 +0.06

Table 2

Values of the migration barriersEX im_→i ,

( ) associated with Ni- and Al-atom jumps to nnn vacancies in NiAl– ab initio calculations[22].

X E_{X α}^{( )}^m_,→_α(eV) E_{X β}^{( )}^m_,_→_β(eV) Average (eV)

Ni 2.76 2.05 2.41

Al 2.42 1.49 1.96

Table 3

Values of the migration-barrier parameters: E_{nn bar}⁺_, ( )A,E_{nnn bar}⁺ _, ( )A, E_{nn bar}⁺_, ( )B andEnnn bar⁺ , ( )B used in the study.

X E_{nn bar}⁺_, ( )X (eV) E_{nnn bar}⁺ _, ( )X(eV)

A 0.4 1.2

B 0.6 0.5

Table 4

Deﬁnitions of normalized (reduced) parameters used in the presentation of the simulation results.

Reduced Parameter Symbol Reference value

Temperature T_red^∗ TO D− ( )S

Temperature Tred T_{O D}⁽^max− ⁾(Fig. 7)

Atom-vacancy pair correlation (C_XV(^μν))red CV(S=0.5,Tred^∗ =1) Tracer diﬀusivity (DX red) D SA( =0.5,Tred=0.47) Activation energy for tracer

diﬀusivity

E D

(A( X))red EA(DA)(S=0.5)

Activation energy for tracer

correlation factor (E_A(f_X⁽^corr⁾))_red EA(DA)(S=0.5) Atomic-jump frequency

w→

( _μ( )^X_ν)red w_α( )^A_→_α(S=0.5,Tred=0.47) Average migration energy 〈E_{X μ}( )^m,_→_ν〉red 〈E_{A α}( )^m,_→_α〉(S=0.5,Tred=0.47)

Fig. 5. Relative chemical potentialsΔμ_{X V}₋ corresponding to the phase equilibria in the A-B-V lattice gas atTred=0.47. Coloured solid circles denote the values determined according to Ref.[29]; open circles represent the values determined by thermodynamic integration[28].

(7)

= −

N_V⁽^αstruct⁾ N_B N_A (20)

= = −

+ + − = −

C N −

N

N N

N N N N

S S

( )

1 2

2(1 )

Vαstruct Vαstruct

B A

B A B A

( ) ( )

(21)

(ii) For structural A-antisites (S > 0.5,CV=0):

− = + =

N N N N N

A Aβstruct 2

B Aβstruct

( ) ( )

(22)

= = −

+ = −

C N

N

N N

N N S

2( )

1

A 2

βstruct Aβstruct

A B

( ) ( )

(23) Fig. 9 shows theS-dependence of the vacancy and antisite concentrations extrapolated toT→0 .K

The curves C_A^{( )}^β ( ,S T→0 )K and C_A⁽^βstruct⁾( )S, as well as

→

C_V^{( )}^α( ,S T 0 )K and C_V⁽^αstruct⁾( )S coincided almost ideally in the range of

>

S 0.5 and 0.45<S<0.5, respectively. This behaviour clearly in- dicated the triple-defect character of the system. In the range of

<

S 0.45 the curveC_V^{( )}^α ( ,S T→0 )K deviated, however, from C_V⁽^αstruct⁾( )S towards lower values ofC_V^{( )}^α which was accompanied by an increase of

the concentration of B-antisites, appearing already atS=0.5(see the inset inFig. 9) and contributing to the compensation of the deficit of A- atoms on the α-sublattice. Remarkably, A-antisites definitely absent belowS=0.5re-appeared atS<0.4. The observed effects, especially the decrease of the vacancy concentration below the value of CV(αstruct)

, indicate that the decrease of A-atom concentration caused a decay of the tendency for TDD of the system.

4.4. Temperature and composition dependence of tracer diﬀusivities of A- and B-atoms

Fig. 10shows the Arrhenius plots of A- and B-atom diﬀusivities evaluated by means of Eq.(13)applied to the KMC-time dependences

〈R_X²( )t〉of the MSD yielded by KMC simulations.

Similarly as in the case ofC_V, CX( )μ

and CXV(μν)

the Arrhenius plots ofDA

andD_Bshowed curvatures close to the order-disorder transition point T_{O D}− , but the linear segments make it possible to evaluate the activation energiesE DA( X)for A- and B-atoms self-diﬀusion (Fig. 11).

Both activation energies showed maximum values close toS=0.5.

WhileE D_A( _A)<E D_A( _B) held forS>0.4, a strong decrease of E D_A( _B) with decreasing S caused that the relationship inverted atS≈0.4. The decrease of E DA( A) for S>0.5 means there is qualitative agreement between the reported simulation results and the corresponding experimental data on Ni-tracer diﬀusion in NiAl[4,5].

Fig. 12 presents the isotherms (DA red) ( )S and (DB red) ( )S corresponding toT_red=0.47 andT_red=0.78. The curves showed the characteristic asymmetric‘V’-shapes with minima located atS=0.43and

=

S 0.5, respectively (see the inset in Fig. 11a). The ‘V’-shape of

D S

( _{B red}) ( ) was deﬁnitely more pronounced and clearly visible in a logarithmic scale (Fig. 11c,d). Besides, the curve increased much stronger with decreasing S than did D( A red) ( )S with increasing S. As a result, the relationshipDA>DBobserved in the range of A-rich binaries inverted atS≈0.4 where the D( _{A red}) ( )S and D( _{B red}) ( )S curves inter- sected. Both features qualitatively corroborate with the experimental results[4,6,7].

Temperature dependences of the positions of the diﬀusivity minima and of the intersection ofD S_A( )andD_B( )S are displayed inFig. 13a and Fig. 6. Sections of the miscibility gap of the A-B-V lattice gas: red solid squares represent the ASB1-Ssystems with an equilibrium vacancy concentration: (a)S=0.31;

(b)S=0.46; (c)S=0.5; (d)S=0.59.

Fig. 7. S-dependence of the reduced‘order-disorder’ transition temperature in the simulated A-B binaries with vacancies.

(8)

b. While the minima of both D S_A( )andD_B( )S shifted towards lower values of S with increasing temperature (Fig. 13a), the location of

=

DA DB remained atS≈0.4 in the whole range of 0.47<Tred<0.8 (Fig. 13b)– which obviously resulted inE D_A( _A)=E D_A( _B)observed at the same value ofS(Fig. 11).

4.5. Elucidation of atomistic origins of the features of A- and B-tracer diﬀusivities

The analysis was based on the atomistic model of self-diﬀusion and

the relationships given by Eqs. (13)–(19) expressing self-diﬀusion coeﬃcients D T SX( , )(observables) in terms of tracer correlation factors f_X⁽^corr⁾( , )T S and atomic-jump frequenciesw_μ^{( )}^X→_ν( , )T S which, in turn, depended on atom-vacancy pair correlations C_XV⁽^μν⁾( , )T S and the average migration barriers 〈E_{X μ}^{( )}^m_, _→_ν〉( , )T S (Eq. (16)). Each one of the above parameters, as well as its composition- and temperature-dependence was independently evaluable by means of MC simulations.

4.5.1. Tracer correlation factors

Fig. 14shows the temperature and composition dependences of the tracer correlation factors f_A⁽^corr⁾and f_B⁽^corr⁾.

The linear parts of the Arrhenius plots of f_X⁽^corr⁾ yielded eﬀective activation energies E fA(X )

corr

( ) traced inFig. 15a against S.Fig. 15b and c show the S-dependence of two relationships between E fA(X )

corr

( ) and the total activation energiesE D_A( _X)for self-diﬀusion (Fig. 11): the diﬀer- ence between both activation energies (Fig. 15b) and the contribution of E f_A(_X⁽^corr⁾)toE D_A( _X)(Fig. 15c).

According to Eq. (14) the activation energy E fA(X )

corr

( ) additively contributes to the activation energyE DA( X)for X-atom tracer-diﬀusion.

The diﬀerenceE D_A( _X)−E f_A(_X⁽^corr⁾)yields, therefore, the part ofE D_A( _X) stemming directly from the kinetics of atomic jumps to vacancies. The graphs inFig. 15c show, in turn, that the contribution of the activation energy E f_A(_X⁽^corr⁾) to the total activation energy E DA( X) for X-tracer diﬀusion never exceeded 30%.

4.5.2. Analysis of D (T, S)X in terms of atomic jump frequencies, atom- vacancy pair-correlations and average migration barriers.

The pure eﬀect of w_μ^{( )}^X→_ν on the diﬀusivities is manifested by the values of D_X( ,S T) evaluated with Eq. (13) divided by the corre- sponding values of f_X⁽^corr⁾( ,S T)evaluated independently with Eq.(15).

Elucidation of the atomistic origin of the observed features of A- and B- atom diﬀusivities follows, in turn, from the analysis of the interrelations Fig.8. Arrhenius plots of (a) total vacancy concentration CV; (b)α-vacancy concentration CV( )α; (c)β-vacancy concentration CV( )βand (d) atom-vacancy pair corre- lations C_XV⁽^μν⁾determined by SGCMC simulations of ASB1-Sbinaries.

Fig. 9. S-dependence of C_V^{( )}^α, C_A^{( )}^β and C_B^{( )}^α extrapolated toT→0K. The solid red line and the dashed black line denote the S-dependences of the concentrations of structural A-antisites and structural A-vacancies, respectively, according to Eqs.(21) and (23).