- Email: [email protected]
MY (α-class)
Solid State Ionics 176 (2005) 2681 – 2688 www.elsevier.com/locate/ssi
Simulation of interdiffusion processes in quasiternary systems of type MX/NX/NY/MY (a-class) Polyzois Kokkonis, Volkmar Leute * Institut fu¨r Physikalische Chemie der Universita¨t Mu¨nster, Corrensstrasse 36, 48149 Mu¨nster, Germany Received 1 April 2005; received in revised form 20 June 2005; accepted 2 September 2005
Abstract Diffusion processes in quasiternary systems are described by a diffusion matrix composed of a matrix covering the kinetic properties and a matrix describing thermodynamic influences by interaction parameters. The interdiffusion processes were simulated by a computer program and the calculated diffusion profiles and reaction paths were compared with earlier reported experimental data. It can be shown that thermodynamics governs the principal reaction paths, whereas the individual diffusion behavior is controlled by the kinetic matrix composed of tracer diffusion coefficients depending on the type of intrinsic doping. Regions of uphill diffusion are determined predominantly by thermodynamic coupling effects. D 2005 Elsevier B.V. All rights reserved. PACS: CdSe; CdTe; HgSe; HgTe Keywords: Interdiffusion; Uphill diffusion; Solid state reactions; Quasiternary systems; Thermodynamics; Strict regular systems; Solid solutions; Simulation processes
1. Introduction Diffusion processes in multicomponent alloys can be designated as solid state reactions. Up to now, in such cases a microscopical description that uses only statistical movements of diffusing particles is practically not possible. Although there are theoretical concepts as master equations or linear response theory the actual problem to calculate reaction paths in thermodynamically real systems including phase boundaries is much too complicated to be solved by these methods. One way to escape this problem is to use the methods of phenomenological thermodynamics that work with fewer variables than statistical dynamics. Moreover, in multicomponent systems coupling effects between diffusion fluxes have to be considered. Thus, a diffusion matrix has to be introduced instead of only one single diffusion coefficient for each component. The most promising way till now to describe the multicomponent diffusion in quasiternary systems is by linear irreversible thermodynamics [1,2]. In this formalism the
* Corresponding author. Tel.: +49 2534 434; fax: +49 251 83 23423. E-mail address: [email protected] (V. Leute). 0167-2738/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.ssi.2005.09.006
diffusion matrix, D, describing the influence of the forces on the fluxes, can be written as a product of two matrices: a kinetic matrix, L, taking the part of the tracer diffusion coefficients, and a thermodynamic matrix, M, playing the part, which the thermodynamic factor plays in binary systems [3]. This matrix, M, is generally described as a function of the interaction parameters that determine the miscibility properties of the solid solution. As will be discussed later, the effect of the offdiagonal components of L is called kinetic coupling and that of the off-diagonal components of M is called thermodynamic coupling. The term Fkinetic_ is justified because the components of L are functions of the tracer diffusion coefficients and the term Fthermodynamic_ because the components of H are the second derivatives of the Gibbs energy. If diffusion processes in quasiternary or higher systems are to be explained, the influence of thermodynamics on diffusion profiles and reaction paths must not be underestimated. The effect of thermodynamics on diffusion processes occurring in quasiternary chalcogenide alloys of type (Mk N1-k )(Xl Y1-l ), called a-class systems, can best be shown by way of solid state reactions of type
MX þ NYYNX þ MY:
ð1Þ
2682
P. Kokkonis, V. Leute / Solid State Ionics 176 (2005) 2681 – 2688
Alloys of the a class are characterized by the fact that the cationic particles M and N at the one side and the anionic particles X and Y at the other side are separated onto two sublattices. As here the cation and anion diffusion fluxes occur on different sublattices, the thermodynamic coupling of fluxes will have more influence on the reaction path than the kinetic coupling. The latter becomes generally more important if different particles use the same diffusion paths. We will investigate in this paper how the observed diffusion behaviour is influenced by thermodynamic properties. To this end we have developed a computer program on the basis of the finite-difference method by which the interdiffusion processes can be simulated and the diffusion profiles and the reaction path can be calculated [4]. These calculated data are compared with experimental results described in previous papers. 2. Theoretical considerations 2.1. The diffusion matrix The driving force for the chemical flux of a component j can be described as: Xj ¼ †lj ¼
X Blj †ck ; Bck k
ð2Þ
where the chemical potential l j of the diffusing component j is given as a function of the component concentrations lj ¼ lj ðc1 ; . . . ; cn Þ:
ð3Þ
Because of the linear relation between fluxes and forces in irreversible thermodynamics [1] X Ji ¼ Lij Xj ð4Þ j
one obtains: Ji ¼
X X j
k
Lij
X Blj †ck ¼ Dik †ck Bck k
ð5Þ
With this equation one has a generalized form of Fick’s first law, in which the terms D ik describe the components of the diffusion coefficient matrix D [5] Dik : ¼
X j
Lij
Blj : Bck
ð6Þ
L ij are the components of the Onsager matrix L describing the kinetic properties of the system, and the derivative fll j / flc k can be interpreted as a component of the Hessian matrix M for the Gibbs energy describing the thermodynamic influence: 2 2 Blj Blj V B g 1 B g 1 ¼V ¼ ¼ ¼ Mjk ð7Þ n Bxj Bxk C Bxj Bxk C Bck Bnk n is the total amount of substance within the considered Volume V and C is the mean molar concentration of the components; x i is the mole fraction of component i. By this
way the diffusion matrix D can be described by the matrix product of the Onsager matrix L and the Hessian matrix M: D¼
1 ðLIMÞ: C
ð8Þ
2.2. Interdiffusion in systems with separated sublattices Diffusion problems for substances with sublattices were ˚ gren [6]. In the present case we made the formulated by A following assumptions: & Only one cationic sublattice and one anionic sublattice are considered; interstitial lattices are not considered as separated sublattices; & the ratio of the numbers of lattice sites on the two sublattices remains constant during diffusion (number fixed reference system); & the particles M and N of the cationic sublattice (c) and X and Y of the anionic sublattice (a) cannot exchange between the two sublattices; & the partial molar volume is the same for all components and independent of composition; the occurrence of interstitial particles does not change the partial molar volume, i.e. the mean molar concentration C is constant. From that we get relations between the fluxes of the i components within a sublattice s: Ai J si = 0. After the reduction to independent fluxes and forces one can write: X c ¼ XMc XNc
and
X a ¼ XXa XYa :
ð9Þ
Schottky [7] and Kro¨ger [8] have pointed out that in thermodynamic descriptions only chemical potentials of building units and not of single defects have to be used and in Ref. [9] we have shown that in quasiternary systems of the a-class, (Mk N1-k )(Xl Y1-l ), the chemical potentials of the difference building units can be expressed by the derivatives of the mean molar Gibbs energy g (k = x M, l = x X): Bg Bg and lðX YÞ ¼ : ð10Þ lðM NÞ ¼ Bk l Bl k If it is assumed that the density of vacancies corresponds everywhere to local equilibrium, and if the flows of particles are taken relative to the fixed parts of the lattice [10] the diffusion equations can be formulated as
B2 g B2 g Bk JM ¼ JN ¼ LMM 2 þ LMX Bk BkBl Bz 2 B g B2 g Bl þ LMX 2 LMM BlBk Bl Bz B2 g B2 g Bk JX ¼ JY ¼ LXM 2 þ LXX Bk BkBl Bz B2 g B2 g Bl þ LXX 2 LXM BkBl Bl Bz
ð11Þ
P. Kokkonis, V. Leute / Solid State Ionics 176 (2005) 2681 – 2688
and from that, using the continuity equation, we obtain a formulation of the 2. Fick’s law with the mole fractions k and l: Bk B Bk Bl ¼ DMM þ DMX ; Bt Bz Bz Bz Bl B Bk Bl ¼ DXM þDXX ; ð12Þ Bt Bz Bz Bz
2.4. The Onsager matrix L The non-diagonal terms: As a result of the second principle of thermodynamics the L matrix can never become negative [11] and consequently LMM 0;
LXX 0;
LMM ILXX
1 ðLMX þ LXM Þ2 : 4
Including the Onsager principle [12] one obtains: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi LMX ¼ LXM ¼ kI LMM ILXX with 1V k V1:
where DAB ¼
2683
1 X B2 g LAK with A; B; KZfM; Kg C K BxK BxB
ð13Þ
are the Fpartial_ diffusion coefficients [9]. This formulation for diffusion in separated sublattices has the advantage that the diffusion processes in each sublattice can be treated independent of each other. In Ref. [9] it was shown that the diffusion processes in the cation sublattice can be interpreted as a interdiffusion of M and N described by an ˜ MN and that in the anion sublattice interdiffusion coefficient D as a interdiffusion of X and Y described by an interdiffusion ˜ XY. coefficient D 2.3. The solution of the diffusion equation by the finite difference method To solve the system of partial diffusion Eq. (12) a discretized grid with (2M + 1) intervals on the distance axis and (N + 1) intervals on the time axis is defined. A ¼ f M ; . . . ; þ M g f0; . . . ; N g:
ð14Þ
The step sizes of the grid are Dz and Dt. Thus, N I Dt corresponds to the total diffusion time t tot. The mole fractions c = k, l to be calculated are defined as a discrete function c on A: cij ¼ cðiDz; jDt Þ j ¼ 0; . . . ; N :
with
i ¼ M; . . . ; þ M
and ð15Þ
As terms for the derivatives, formulated by finite differences, we use: ci;jþ1 ci;j Bc ’ Bt Dt
ð16Þ
and Di;j þ Diþ1;j ciþ1;j ci;j Di;j þ Di1;j ci;j ci1;j B Bc D ’ : Bz Bz 2ð DzÞ2
ð17Þ
For each term of the sums in Eq. (12) a term, as described by Eq. (17), has to be used.
ð18Þ
ð19Þ
The parameter k describes ternary effects caused by kinetic coupling and depends generally on composition and temperature. But if, as we assumed, the two independent fluxes J M and J X (Eq. (11) occur in two separate sublattices, the mobility of defects of the one sublattice should not influence the mobility of the diffusing particles in the other sublattice and vice versa. That is, a kinetic coupling should not occur (L MX = L XM = 0 Y k = 0) provided the defect density is low enough. At high defect densities, however, the potential surface around the defects could influence the potential barriers for the diffusing particles in the other sublattice and then a kinetic coupling should be observable. As the defect densities of the binary components of the given system are rather small, we take k = 0 for the following calculations. The diagonal terms: According to Darken’s equation for binary thermodynamically ideal systems [13] D˜ MN;ideal ¼ kDN þ ð1 k ÞDM
ð20Þ
we use a model, in which we assume that the ideal interdiffusion coefficients in all quasibinary subsystems, e.g. in (MX)k / (NX)1-k , can be described as being linearly dependent on composition. This means that the diffusion coefficients in Eq. (20) become equal to the extrapolated tracer diffusion coefficients D N(MX) of N in pure MX or D M(NX) of M in pure NX. ˜ AB,ideal = RT[(x B / C A)L A + In analogy to the equation D (x B / C B)L B] for binary systems, as explained in [14], we use the ˜ MN,ideal = RT I L MM / C and D ˜ XY,ideal = similar relations D RT I L XX / C for the calculation of the diagonal coefficients of the Onsager matrix L in our quasiternary system: LMM 1 fl ½kIDN ðMXÞ þ ð1 k ÞIDM ðNXÞ ¼ RT C þ ð1 l Þ½kIDN ðMYÞ þ ð1 k ÞIDM ðNYÞ g
ð21Þ
and LXX 1 fk ½lIDY ðMXÞ þ ð1 l ÞIDX ðMYÞ ¼ RT C þ ð1 k Þ½lIDY ðNXÞ þ ð1 l ÞIDX ðNYÞ g:
ð22Þ
The extrapolated diffusion coefficients can experimentally be determined by measuring the composition dependent ˜ along a quasibinary edge of the interdiffusion coefficient D quasiternary system and by subsequent extrapolation of these values up to the pure binary components.
2684
P. Kokkonis, V. Leute / Solid State Ionics 176 (2005) 2681 – 2688
Although the relations between the interdiffusion coefficient for a sublattice and the corresponding diagonal Onsager coefficient cannot exactly be deduced from theoretical reasons, they yield plausible procedures that will allow us to estimate composition dependent diagonal Onsager coefficients needed for the calculation of reaction paths in good accordance with the experimental results.
E ðk; lÞ M12
2.5. The Hessian matrix M To calculate the components of the Hessian matrix we need the mean molar Gibbs energy g (k, l). Analogously to the splitting of g also the M matrix can be split into a standard term, an ideal mixing term and an excess term: g ¼ g 0 þ g m þ g E Y M ¼ M 0 þ Mm þ ME : ð23Þ The standard term M0 can be omitted because all second derivatives of g 0 vanish. The ideal mixing term g m is restricted to the ideal mixing entropy as the ideal mixing enthalpy is zero. h i gm ¼ RT Iln k k ð1 k Þð1k Þ l l ð1 l Þð1lÞ : ð24Þ Thus m ¼ RT I M11
1 1 m m m ; M22 ¼ RT I ; M12 ¼ M21 ¼ 0: k k2 l l2 ð25Þ
The term for the excess energy of a quasiternary a-class system (Mk N1-k )(Xl Y1-l ) was derived in [15]: gE ðk; l Þ ¼ k ð1 k Þ½aðMNÞ þ kIbðMNÞ þ lIcðMNÞ þ kIlIdðMNÞ þ l ð1 l Þ½aðXYÞ þ lIbðXYÞ þ kIcðXYÞ þ lIkIdðXYÞ þ k ð1 l ÞIDR G0 :
ð26Þ
As also shown in [15], the quasiternary interaction parameters a(i), h(i), g(i) and y(i) in Eq. (26) can be calculated from the interaction parameters determined for the four quasibinary subsystems. Treating these binary interaction parameters as depending on temperature means that an excess entropy is considered. The ternary parameter D R G 0 =D R H 0
(a)
T I D R S 0 is the molar standard Gibbs energy for the displacement reaction (1). The components of the excess term ME of the Hessian matrix describing the deviation from ideal behaviour read as follows: E M11 ðk; l Þ ¼ ð2 6k ÞðbðMNÞ þ lIdðMNÞÞ 2ðaðMNÞ þ lIcðMNÞÞ E ðk; l Þ ¼ ð1 2k ÞIcðMNÞ þ 2k 3k 2 IdðMNÞ ¼ M21 þ ð1 2l ÞIcðXYÞ þ 2l 3l2 IdðXYÞ DR G0 E M22 ðk; lÞ ¼ ð2 6l ÞðbðXYÞ þ kIdðXYÞÞ 2ðaðXYÞ þ kIcðXYÞÞ
ð27Þ The ME matrix is symmetrical because of M12E = M21E. The non-diagonal term, being important for the ternary effects, depends on D R G 0 and its dependence on composition is only determined by the parameters c and d. Therefore, in symmetrical systems, where neither the interaction between M and N depends on l (c(MN) = 0), d(MN) = 0) nor the interaction between X and Y on k (c(XY) = 0), d(XY) = 0), the whole non-diagonal term will be reduced to M12 = M21 = D R G 0 [15]. That means that in these symmetrical cases D R G 0 is the only parameter that can cause a thermodynamical coupling between the components of the two sublattices. 3. Simulations of interdiffusion in strict regular quasiternary systems In strict regular systems there is per definition no temperature or composition dependence of the interaction parameters. As shown in Ref. [16] the molar standard Gibbs energy of the reaction MX + NY Y NX + MY (1) is identical with the interaction parameter a(MY), i.e. D R G 0 = a(MY). Thus, the excess part of the M matrix can here be written as: 2aðMNÞ aðMNÞ E : ð28Þ M ¼ aðMYÞ 2aðXYÞ The influence of the thermodynamic properties on the diffusion behaviour of quasiternary systems can best be
(b)
Fig. 1. Reaction paths with a(MY) > 0 for D 2 = 10 I D 1 (a) and D 1 = 10 I D 2 (b).
P. Kokkonis, V. Leute / Solid State Ionics 176 (2005) 2681 – 2688 Table 1 Interaction parameters Quasibinary system i
(Hgk Cd1k )Te
(Hgk Cd1k )Se
Hg(Tel Se1l )
Cd(Tel Se1l )
a(i) / (J/mol)
5900
5900
6300
6300
D R(1)G 0 = a(MY) = 20 500
J/mol;
D R(1)S 0 å 0.
realized by the so-called reaction path C(k(z), l(z)) that can be constructed from the mole fraction profiles k(z) and l(z) [9]. Generally, to explain the diffusion behaviour one uses the projection of the reaction path onto an iso-thermic section l(k) of the phase diagram of the quasiternary system (cf. Fig. 1). The most important feature that determines the shape of a reaction path is a spinodal miscibility gap, because each reaction path touching the spinodal curve crosses the region of instability along a tie line. For a strict regular quasiternary system the isothermal spinodal curve is given by the equation xy 2aðXYÞIx þ 2aðMNÞIy þ 4aðMNÞIaðXYÞ aðMYÞ2 ¼ 0 ð29Þ with x := RT / [k(1 k)] and y := RT / [l(1 l)]. In Ref. [16] it was shown that this means that the shape of the spinodal curve at a given temperature only depends on the magnitude of a(MY), independent of whether the displacement reaction MX + NY Y NX + MY is exo- or endothermic. The shape of the equilibrium curve, however, depends on the sign of a(MY). To show that thermodynamics influences diffusion even in this most simple case, we make very simple assumptions for the L matrix. The kinetic coupling effects shall be negligible which means that the L matrix must be diagonal, i.e. k = 0. Furthermore, the L matrix is taken as independent of composition (k, l) which means that all four diffusion coefficients in Eq. (21) or Eq. (22) must have identical values, i.e. there are only two diffusion coefficients one for the cationic (D 1) and the other for the anionic components (D 2). Fig. 1a and b shows the reaction paths for an endothermic reaction MX + NY Y NX + MY and its reverse exothermic reaction NX + MY Y MX + NY on condition that D 2 = 10 I D 1
2685
(a) or D 1 = 10 I D 2 (b). These reaction paths were simulated by computer calculations using the values for the quasibinary interaction parameters as shown in Table 1. In [4] we have shown that by changing the sign of a(MY) the reaction path for a strict regular system as well as the tie lines and the spinodal curve are reflected at the axis k = 0.5. Therefore only the case a(MY) > 0 is treated here explicitly. Within the spinodal curve the reaction path runs always along the tie lines. Therefore, independent of the diffusion coefficients, the reaction path for an endothermic reaction runs always nearly along the diagonal of the phase square. An exothermic reaction, however, shows quite another behaviour. Here, the reaction path also crosses the spinodal curve via a tie line, but it starts more or less along quasibinary edges until a tie line is reached. If the diffusion coefficient for the chalcogens X and Y is higher than that for the constituents of the metal sublattice M and N (Fig. 1a) the reaction path runs predominantly along the chalcogen edges M(X,Y) and N(X,Y), in the reverse case (Fig. 1b) along the edges (M,N)X and (M,N)Y. Whereas the jump across the miscibility gap is forced by thermodynamics, the shape of the reaction path within the homogeneous region is determined predominantly by kinetics, i.e. by the ratio of the interdiffusion coefficients for the two sublattices. This kinetic effect becomes also apparent near the pure binary components of the endothermic reaction. In [17] we had already demonstrated that a Z-shaped reaction path, as shown in Fig. 1b, is retained even at such high temperatures, where the spinodal miscibility gap does no longer exist. At temperatures above the critical temperature the thermodynamic interaction causes a region of uphill diffusion caused by the off-diagonal part of the Hessian matrix. This behaviour is quite analogous to the influence of the thermodynamic factor on interdiffusion at temperatures above spinodal miscibility gaps in quasibinary systems [3]. Till now the initial conditions (t = 0) for the simulation of the displacement reactions have been the following: ðkV ¼ 0; lV ¼ 1Þ; ðkþV ¼ 1; lþV ¼ 0Þ or ðkV ¼ 1; lV ¼ 0Þ; ðkþV ¼ 0; lþV ¼ 1Þ
Fig. 2. Reaction paths with D 1 = 10 I D2 for a(MY) > 0.
2686
P. Kokkonis, V. Leute / Solid State Ionics 176 (2005) 2681 – 2688
(a)
(c)
(b)
(d)
Fig. 3. Diffusion profiles k(z) and l(z) corresponding to reaction paths selected from Fig. 2 (a(MY) > 0).
But the Hessian matrix M is not influenced by the initial conditions, thus it remains unchanged also for so-called quasibinary initial conditions: ðkV ¼ c; lV ¼ 1Þ; ðkþV ¼ c; lþV ¼ 0Þ; ðFig: 2aÞ or ðkV ¼ 0; lV ¼ cÞ; ðkþV ¼ 1; lþV ¼ cÞ; ðFig: 2bÞ
with
0 V c V 1: The simulations of reaction paths with quasibinary initial conditions, as shown in Fig. 2a and b, yield pronounced deviations from quasibinary sections. These deviations are caused by the ternary interaction parameter a(MY). As long as the reaction path does not intersect the spinodal miscibility gap one obtains reaction paths corresponding to diffusion profiles with regions of uphill diffusion for those components with identical boundary conditions on both ends of the profile; e.g. for interdiffusion of Hg and Cd in Fig. 3a or of Te and Se in Fig. 3c. If the reaction path crosses the spinodal miscibility gap, the diffusion profiles show reverse concentration steps at the phase boundary as for example Hg in Fig. 3b. The comparison of Fig. 3a and b shows that, if the phase diagram were not known, one could hardly distinguish between continuous profiles with uphill diffusion and discontinuous profiles with inverse concentration steps at phase boundaries. Even more complicated are situations where the reaction path cannot cross the spinodal miscibility gap along a tie line because it reaches the spinodal curve near one of the critical points for spinodal demixing as in Fig. 2b for l = 0.75. In such cases one obtains an oscillating behaviour until the reaction path can leave the unstable region (Fig. 3d). 4. The system (Hgk Cd1-k )(Tel Se1-l ) 4.1. Results of computer simulations To show that the kinetic influence on solid state reactions can be reproduced correctly by the computer simulations, diffusion profiles for the displacement reaction HgTe + CdSe Y CdTe + HgSe (1) were simulated with different sets of
Table 2 Experimentally determined interdiffusion coefficients at 723 K [19]; (k, l) are the compositions at which the interdiffusion coefficients were measured, the subscripts Fme_ and Fch_ mean metal or chalcogen saturated ˜ meexp / (cm2/s) ˜ chexp / (cm2/s) Quasibinary subsystem (k, l) D D (Hgk Cd1k )Te (Hgk Cd1k )Te (Hgk Cd1k )Se (Hgk Cd1k )Se Hg(Tel Se1l ) Hg(Tel Se1l ) Cd(Tel Se1l ) Cd(Tel Se1l )
(0.1, (0.9, (0.9, (0.1, (1.0, (1.0, (0.0, (0.0,
1.0) 1.0) 0.0) 0.0) 0.9) 0.1) 0.1) 0.9)
2 I 10 14 3 I 10 12 2 I 10 10 6 I 10 13 3 I 10 14 1.5 I 10 14 2 I 10 17 5 I 10 16
1.4 I 10 13 3 I 10 12 2 I 10 11 3 I 10 13 4 I 10 10 2 I 10 10 – 1 I10 13
diffusion coefficients. The same reaction was already experimentally investigated in dependence on intrinsic doping in [18,19]. The interaction parameters for this system as derived in Ref. [18] are listed in Table 1. We tried to reproduce the diffusion profiles and reaction paths for metal saturated and chalcogen saturated samples at 723 K and for samples with the composition corresponding to minimum vapor pressure at 673 K. To find suitable sets of tracer diffusion coefficients, as defined by Eq. (21) or (22), we used as starting values the interdiffusion coefficients for compositions near the corners of the phase diagram, which were experimentally determined in Ref. [19] for reactions between metal saturated as well as between chalcogen saturated samples (Table 2). The tracer diffusion coefficients estimated by adjusting the calculated diffusion profiles and reaction paths to the experimentally determined data are listed in Table 3. The calculated reaction paths and diffusion profiles are plotted in Figs. 4a –c and 5a –c; they correspond rather good with the experimentally determined ones as shown in [18,19]. It is astonishing that the rather simple assumptions, used for modeling the Onsager matrix, allow to reproduce the experimentally determined diffusion profiles and reaction paths in all essential details. For the simple model used here, all kinetic and thermodynamic data that are needed for the simulation of the diffusion profiles and reaction paths can be taken from experiments on the quasibinary subsystems: & The molar standard Gibbs energy for the displacement reaction D R G 0 follows principally from the standard chemical potentials of the pure binary components. & The ternary interaction parameters (cf. Eq. (26)) can be calculated from the interaction parameters of the quasibinary Table 3 Tracer diffusion coefficients of minor components in the binary compounds determined by simulation; the subscripts Fme_, Fch_ and Fmin_ mean metal saturated, chalcogen saturated or minimal vapor pressure ˜ me / (cm2/s) ˜ ch / (cm2/s) ˜ min / (cm2/s) Tracer diffusion coeficient D D D Hg in CdTe Cd in HgTe Cd in HgSe Hg in CdSe Se in HgTe Te in HgSe Te in CdSe Se in CdTe
1 I10 13 1 I10 12 4 I 10 11 1 I10 12 1 I10 13 1.5 I 10 14 2 I 10 17 1 I10 15
5 I 10 12 3 I 10 12 1 I10 12 3 I 10 13 4 I 10 10 2 I 10 10 1 I10 12 1 I10 12
1 I10 12 1 I10 13 1 I10 13 1 I10 13 5 I 10 12 1 I10 12 1 I10 12 1 I10 12
P. Kokkonis, V. Leute / Solid State Ionics 176 (2005) 2681 – 2688
(a)
2687
(b)
(c)
Fig. 4. Reaction paths for (a) chalcogen and (b) metal saturation and for (c) minimum vapor pressure. The dashed area describes the spinodal miscibility gap with tie lines and the inner curve is the spinodal curve.
subsystems even for nonregular systems. The spinodal miscibility gap– being in accordance with the experimental phase diagram – was calculated with the same set of interaction parameters. & The diffusion coefficients needed for constructing the Onsager matrix according to Eq. (21) or (22) correspond to the extrapolated interdiffusion coefficients of the quasibinary subsystems (Table 2). According to Darken’s equation one obtains by extrapolation of the pure binary components the tracer diffusion coefficients of the minor component, e.g. for the quasibinary subsystem (M,N)X the diffusion coefficients for N in MX and M in NX.
& The coefficient k (Eq. (19)), responsible for kinetically caused ternary diffusion effects, can in most systems with separated sublattices be taken as zero. Corresponding to earlier reported experimental findings [19] the computer simulations yield continuous reaction paths for chalcogen (Fig. 4a) as well as for metal saturated (Fig. 4b) reaction couples. For a reaction couple with the composition of minimal vapor pressure the simulations yield – also corresponding to the experiments [18] – a reaction path crossing the spinodal miscibility gap. 4.2. Uphill diffusion
(a)
(b)
(c)
The diffusion profiles belonging to the continuous reaction paths show distinct zones of uphill diffusion. Starting at the compound HgSe the reaction path runs at first in that direction, where the interdiffusion coefficients are the highest, i.e. for the chalcogen saturated case along the edge Hg(Se,Te) (Fig. 4a) or for the metal saturated case along the edge (Hg, Cd)Se (Fig. 4b). In the metal saturated case the decrease of the Hg concentration is stopped, when the chalcogen exchange starts. In the same measure as the chalcogen exchange increases, the uphill diffusion of Hg gets promoted (cf. Fig. 5b). The essential term for this uphill diffusion is the partial diffusion coefficient D MX of Eq. (12) for M = Hg and X = Te: DHgTe
Fig. 5. Diffusion profiles k(z) and l(z) for (a) chalcogen (t = 1 h) and (b) metal saturation (t = 1 h) and for (c) minimum vapor pressure (t = 5 h).
1 B2 g B2 g LHgHg ¼ þ LHgTe C BxHg BxTe BxT e BxTe ¼
1 LHgHg I aHgTe : C
ð30Þ
2688
P. Kokkonis, V. Leute / Solid State Ionics 176 (2005) 2681 – 2688
For the given system this last equation becomes so simple as L HgTe = 0 because of k = 0 and as in strictly regular systems the off-diagonal terms of the matrix (Eq. (28)) consist only of the ternary interaction parameter a HgTe. Obviously, ternary effects like uphill diffusion can occur without kinetic coupling effects. In most cases uphill diffusion can be explained as caused only by thermodynamics, here for example by a HgTe. The partial diffusion coefficient D HgTe in Eq. (30) has a negative sign because of the positive sign of a HgTe (Table 1). Therefore the Hg flux according to Eq. (31) and Fig. 5b can remain negative even if the gradient flk / flz for this flux becomes negative, provided the slope of the Te profile flk / flz is high enough and also negative. Bk Bl Bk þ DHgTe : ð31Þ JHg ¼ C DHgHg ¼ C D˜ M Bz Bz Bz The sign of the ternary interdiffusion coefficient for the metal sublattice [9], defined by D˜ M : ¼ DHgHg þ DHgTe Bl=Bk;
ð32Þ
depends on the slope of the reaction path, fll / flk. An uphill diffusion of the cationic components demands (Eq. (31)) a ˜ M and thus, because of D HgHg > 0 (Eq. negative value for D (18)), it demands opposite signs for DHgTe and fll / flk. Furthermore, the magnitude of the second term in Eq. (32) must be high enough to over-compensate the first term. This can be easier attained the higher the slope of the reaction path. The discussion of uphill diffusion for the chalcogen saturated case has to proceed analogously. 5. Conclusions The calculations of interdiffusion profiles and reaction paths for solid state reactions in a-class systems have demonstrated that the experimentally observed phenomena can be explained by a rather simple model. Whether the reaction path has to cross miscibility gaps or runs along continuous solid solutions depends on the Hessian matrix. Especially ternary effects like uphill diffusion can be understood as caused predominantly by thermodynamic coupling effects initiated by the interaction parameter a(MY) being identical with D R G 0 of the displacement reaction MX + NY Y NX + MY. The influences of intrinsic doping on diffusion are described by the Onsager matrix. In those calculations where the reaction path crosses the spinodal miscibility gap a problem arises, because the region between the spinodal curve and the equilibrium curve of the spinodal miscibility gap is still stable with regard to diffusion [20]. Therefore, a reaction path that is calculated with the described simulation model does not become unstable before it reaches the spinodal curve. Under conditions of local equilibrium, however, the reaction path should end already at the equilibrium curve. To solve this problem the precipitation processes within the metastable region would have to be considered, i.e. nucleation kinetics would have to be
included into the diffusion algorithms, a nearly unsolvable task. But compared with experiments, this problem is not so serious at a glance. A diffusion profile crossing a phase boundary, determined experimentally by an electron microprobe, can hardly be distinguished from a calculated diffusion profile, because of the sharp concentration changes near the interface. Calculated profiles yield in this small region only a few calculated points and measurements with a microprobe do also not yield sharp concentration steps, but a more or less continuous course, because of the restricted lateral resolution of the measured X-ray signals. Thus, considering that one cannot even be sure that during diffusion processes local equilibrium at the interface is really established, the calculated profiles correspond unexpectedly well with the experimental results. All things considered, the simulations of reaction paths and concentration profiles, on the basis of the described model, allow us to understand ternary effects like uphill diffusion and to estimate sets of diffusion coefficients by which the essentials of measured concentration profiles can be reproduced. (The computer programs for the simulation processes can be downloaded from the Internet [21].) References [1] P. Glansdorff, I. Prigogine, Thermodynamic Theory of Structure, Stability and Fluctuations, Wiley-Interscience, London, 1971, p. 30. [2] A.R. Allnatt, A.B. Lidiard, Atomic Transport in Solids, Cambridge Univ. Press, Cambridge, 1993. [3] V. Leute, Solid State Ion. 164 (2003) 159. [4] A. Polyzois, Kokkonis, Die Rolle der Thermodynamik fu¨r die Deutung von Kopplungseffekten bei mehrkomponentigen Diffusionsprozessen Thesis (1997) Mu¨nster. [5] J.S. Kirkaldy, D.J. Young, Diffusion in the Condensed State, The Institute of Metals, London, p. 154f. ˚ gren, J. Phys. Chem. Solids 43 (1982) 421. [6] J. A [7] W. Schottky, in: W. Schottky (Ed.), Halbleiterprobleme, vol. 4, Vieweg, Braunschweig, 1958, p. 235. [8] F.A. Kro¨ger, The Chemistry of Imperfect Crystals, Vol. 1-3, North Holland, Amsterdam, 1974. [9] V. Leute, Solid State Ion. 17 (1985) 185. [10] A.R. Allnatt, A.B. Lidiard, Atomic Transport in Solids, Cambridge Univ. Press, Cambridge, 1993, p. 188. [11] S.R. de Groot, P. Mazur, Grundlagen der Thermodynamik Irreversibler Prozesse, Bibliograph. Inst., Mannheim, 1969, p. 35. [12] L. Onsager, Phys. Rev. 37 (1931) 405. [13] A.R. Allnatt, A.B. Lidiard, Atomic Transport in Solids, Cambridge Univ. Press, Cambridge, 1993, p. 197. [14] J.S. Kirkaldy, D.J. Young, Diffusion in the Condensed State, The Inst. of Metals, London, 1987, p. 142. [15] V. Leute, Ber. Bunsenges. Phys. Chem. 93 (1989) 7. [16] V. Leute, Ber. Bunsenges. Phys. Chem. 80 (1976) 866. [17] P.A. Kokkonis, V. Leute, Defect Diffus. Forum 143/147 (1997) 1159. [18] V. Leute, H.M. Schmidtke, W. Stratmann, Ber. Bunsenges. Phys. Chem. 86 (1982) 732. [19] V. Leute, H.M. Schmidtke, W. Stratmann, A. Kalb, Ber. Bunsenges. Phys. Chem. 87 (1983) 483. [20] I. Prigogine, R. Defay, Chemische Thermodynamik, VEB Deutscher Verlag fu¨r Grundstoffindustrie, Leipzig, 1962, p. 243. [21] http://www.uni-muenster.de/Chemie.pc/Leute/Welcome.html.
Simulation of interdiffusion processes in quasiternary systems of type MX/NX/NY/MY (a-class) Polyzois Kokkonis, Volkmar Leute * Institut fu¨r Physikalische Chemie der Universita¨t Mu¨nster, Corrensstrasse 36, 48149 Mu¨nster, Germany Received 1 April 2005; received in revised form 20 June 2005; accepted 2 September 2005
Abstract Diffusion processes in quasiternary systems are described by a diffusion matrix composed of a matrix covering the kinetic properties and a matrix describing thermodynamic influences by interaction parameters. The interdiffusion processes were simulated by a computer program and the calculated diffusion profiles and reaction paths were compared with earlier reported experimental data. It can be shown that thermodynamics governs the principal reaction paths, whereas the individual diffusion behavior is controlled by the kinetic matrix composed of tracer diffusion coefficients depending on the type of intrinsic doping. Regions of uphill diffusion are determined predominantly by thermodynamic coupling effects. D 2005 Elsevier B.V. All rights reserved. PACS: CdSe; CdTe; HgSe; HgTe Keywords: Interdiffusion; Uphill diffusion; Solid state reactions; Quasiternary systems; Thermodynamics; Strict regular systems; Solid solutions; Simulation processes
1. Introduction Diffusion processes in multicomponent alloys can be designated as solid state reactions. Up to now, in such cases a microscopical description that uses only statistical movements of diffusing particles is practically not possible. Although there are theoretical concepts as master equations or linear response theory the actual problem to calculate reaction paths in thermodynamically real systems including phase boundaries is much too complicated to be solved by these methods. One way to escape this problem is to use the methods of phenomenological thermodynamics that work with fewer variables than statistical dynamics. Moreover, in multicomponent systems coupling effects between diffusion fluxes have to be considered. Thus, a diffusion matrix has to be introduced instead of only one single diffusion coefficient for each component. The most promising way till now to describe the multicomponent diffusion in quasiternary systems is by linear irreversible thermodynamics [1,2]. In this formalism the
* Corresponding author. Tel.: +49 2534 434; fax: +49 251 83 23423. E-mail address: [email protected] (V. Leute). 0167-2738/$ - see front matter D 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.ssi.2005.09.006
diffusion matrix, D, describing the influence of the forces on the fluxes, can be written as a product of two matrices: a kinetic matrix, L, taking the part of the tracer diffusion coefficients, and a thermodynamic matrix, M, playing the part, which the thermodynamic factor plays in binary systems [3]. This matrix, M, is generally described as a function of the interaction parameters that determine the miscibility properties of the solid solution. As will be discussed later, the effect of the offdiagonal components of L is called kinetic coupling and that of the off-diagonal components of M is called thermodynamic coupling. The term Fkinetic_ is justified because the components of L are functions of the tracer diffusion coefficients and the term Fthermodynamic_ because the components of H are the second derivatives of the Gibbs energy. If diffusion processes in quasiternary or higher systems are to be explained, the influence of thermodynamics on diffusion profiles and reaction paths must not be underestimated. The effect of thermodynamics on diffusion processes occurring in quasiternary chalcogenide alloys of type (Mk N1-k )(Xl Y1-l ), called a-class systems, can best be shown by way of solid state reactions of type
MX þ NYYNX þ MY:
ð1Þ
2682
P. Kokkonis, V. Leute / Solid State Ionics 176 (2005) 2681 – 2688
Alloys of the a class are characterized by the fact that the cationic particles M and N at the one side and the anionic particles X and Y at the other side are separated onto two sublattices. As here the cation and anion diffusion fluxes occur on different sublattices, the thermodynamic coupling of fluxes will have more influence on the reaction path than the kinetic coupling. The latter becomes generally more important if different particles use the same diffusion paths. We will investigate in this paper how the observed diffusion behaviour is influenced by thermodynamic properties. To this end we have developed a computer program on the basis of the finite-difference method by which the interdiffusion processes can be simulated and the diffusion profiles and the reaction path can be calculated [4]. These calculated data are compared with experimental results described in previous papers. 2. Theoretical considerations 2.1. The diffusion matrix The driving force for the chemical flux of a component j can be described as: Xj ¼ †lj ¼
X Blj †ck ; Bck k
ð2Þ
where the chemical potential l j of the diffusing component j is given as a function of the component concentrations lj ¼ lj ðc1 ; . . . ; cn Þ:
ð3Þ
Because of the linear relation between fluxes and forces in irreversible thermodynamics [1] X Ji ¼ Lij Xj ð4Þ j
one obtains: Ji ¼
X X j
k
Lij
X Blj †ck ¼ Dik †ck Bck k
ð5Þ
With this equation one has a generalized form of Fick’s first law, in which the terms D ik describe the components of the diffusion coefficient matrix D [5] Dik : ¼
X j
Lij
Blj : Bck
ð6Þ
L ij are the components of the Onsager matrix L describing the kinetic properties of the system, and the derivative fll j / flc k can be interpreted as a component of the Hessian matrix M for the Gibbs energy describing the thermodynamic influence: 2 2 Blj Blj V B g 1 B g 1 ¼V ¼ ¼ ¼ Mjk ð7Þ n Bxj Bxk C Bxj Bxk C Bck Bnk n is the total amount of substance within the considered Volume V and C is the mean molar concentration of the components; x i is the mole fraction of component i. By this
way the diffusion matrix D can be described by the matrix product of the Onsager matrix L and the Hessian matrix M: D¼
1 ðLIMÞ: C
ð8Þ
2.2. Interdiffusion in systems with separated sublattices Diffusion problems for substances with sublattices were ˚ gren [6]. In the present case we made the formulated by A following assumptions: & Only one cationic sublattice and one anionic sublattice are considered; interstitial lattices are not considered as separated sublattices; & the ratio of the numbers of lattice sites on the two sublattices remains constant during diffusion (number fixed reference system); & the particles M and N of the cationic sublattice (c) and X and Y of the anionic sublattice (a) cannot exchange between the two sublattices; & the partial molar volume is the same for all components and independent of composition; the occurrence of interstitial particles does not change the partial molar volume, i.e. the mean molar concentration C is constant. From that we get relations between the fluxes of the i components within a sublattice s: Ai J si = 0. After the reduction to independent fluxes and forces one can write: X c ¼ XMc XNc
and
X a ¼ XXa XYa :
ð9Þ
Schottky [7] and Kro¨ger [8] have pointed out that in thermodynamic descriptions only chemical potentials of building units and not of single defects have to be used and in Ref. [9] we have shown that in quasiternary systems of the a-class, (Mk N1-k )(Xl Y1-l ), the chemical potentials of the difference building units can be expressed by the derivatives of the mean molar Gibbs energy g (k = x M, l = x X): Bg Bg and lðX YÞ ¼ : ð10Þ lðM NÞ ¼ Bk l Bl k If it is assumed that the density of vacancies corresponds everywhere to local equilibrium, and if the flows of particles are taken relative to the fixed parts of the lattice [10] the diffusion equations can be formulated as
B2 g B2 g Bk JM ¼ JN ¼ LMM 2 þ LMX Bk BkBl Bz 2 B g B2 g Bl þ LMX 2 LMM BlBk Bl Bz B2 g B2 g Bk JX ¼ JY ¼ LXM 2 þ LXX Bk BkBl Bz B2 g B2 g Bl þ LXX 2 LXM BkBl Bl Bz
ð11Þ
P. Kokkonis, V. Leute / Solid State Ionics 176 (2005) 2681 – 2688
and from that, using the continuity equation, we obtain a formulation of the 2. Fick’s law with the mole fractions k and l: Bk B Bk Bl ¼ DMM þ DMX ; Bt Bz Bz Bz Bl B Bk Bl ¼ DXM þDXX ; ð12Þ Bt Bz Bz Bz
2.4. The Onsager matrix L The non-diagonal terms: As a result of the second principle of thermodynamics the L matrix can never become negative [11] and consequently LMM 0;
LXX 0;
LMM ILXX
1 ðLMX þ LXM Þ2 : 4
Including the Onsager principle [12] one obtains: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi LMX ¼ LXM ¼ kI LMM ILXX with 1V k V1:
where DAB ¼
2683
1 X B2 g LAK with A; B; KZfM; Kg C K BxK BxB
ð13Þ
are the Fpartial_ diffusion coefficients [9]. This formulation for diffusion in separated sublattices has the advantage that the diffusion processes in each sublattice can be treated independent of each other. In Ref. [9] it was shown that the diffusion processes in the cation sublattice can be interpreted as a interdiffusion of M and N described by an ˜ MN and that in the anion sublattice interdiffusion coefficient D as a interdiffusion of X and Y described by an interdiffusion ˜ XY. coefficient D 2.3. The solution of the diffusion equation by the finite difference method To solve the system of partial diffusion Eq. (12) a discretized grid with (2M + 1) intervals on the distance axis and (N + 1) intervals on the time axis is defined. A ¼ f M ; . . . ; þ M g f0; . . . ; N g:
ð14Þ
The step sizes of the grid are Dz and Dt. Thus, N I Dt corresponds to the total diffusion time t tot. The mole fractions c = k, l to be calculated are defined as a discrete function c on A: cij ¼ cðiDz; jDt Þ j ¼ 0; . . . ; N :
with
i ¼ M; . . . ; þ M
and ð15Þ
As terms for the derivatives, formulated by finite differences, we use: ci;jþ1 ci;j Bc ’ Bt Dt
ð16Þ
and Di;j þ Diþ1;j ciþ1;j ci;j Di;j þ Di1;j ci;j ci1;j B Bc D ’ : Bz Bz 2ð DzÞ2
ð17Þ
For each term of the sums in Eq. (12) a term, as described by Eq. (17), has to be used.
ð18Þ
ð19Þ
The parameter k describes ternary effects caused by kinetic coupling and depends generally on composition and temperature. But if, as we assumed, the two independent fluxes J M and J X (Eq. (11) occur in two separate sublattices, the mobility of defects of the one sublattice should not influence the mobility of the diffusing particles in the other sublattice and vice versa. That is, a kinetic coupling should not occur (L MX = L XM = 0 Y k = 0) provided the defect density is low enough. At high defect densities, however, the potential surface around the defects could influence the potential barriers for the diffusing particles in the other sublattice and then a kinetic coupling should be observable. As the defect densities of the binary components of the given system are rather small, we take k = 0 for the following calculations. The diagonal terms: According to Darken’s equation for binary thermodynamically ideal systems [13] D˜ MN;ideal ¼ kDN þ ð1 k ÞDM
ð20Þ
we use a model, in which we assume that the ideal interdiffusion coefficients in all quasibinary subsystems, e.g. in (MX)k / (NX)1-k , can be described as being linearly dependent on composition. This means that the diffusion coefficients in Eq. (20) become equal to the extrapolated tracer diffusion coefficients D N(MX) of N in pure MX or D M(NX) of M in pure NX. ˜ AB,ideal = RT[(x B / C A)L A + In analogy to the equation D (x B / C B)L B] for binary systems, as explained in [14], we use the ˜ MN,ideal = RT I L MM / C and D ˜ XY,ideal = similar relations D RT I L XX / C for the calculation of the diagonal coefficients of the Onsager matrix L in our quasiternary system: LMM 1 fl ½kIDN ðMXÞ þ ð1 k ÞIDM ðNXÞ ¼ RT C þ ð1 l Þ½kIDN ðMYÞ þ ð1 k ÞIDM ðNYÞ g
ð21Þ
and LXX 1 fk ½lIDY ðMXÞ þ ð1 l ÞIDX ðMYÞ ¼ RT C þ ð1 k Þ½lIDY ðNXÞ þ ð1 l ÞIDX ðNYÞ g:
ð22Þ
The extrapolated diffusion coefficients can experimentally be determined by measuring the composition dependent ˜ along a quasibinary edge of the interdiffusion coefficient D quasiternary system and by subsequent extrapolation of these values up to the pure binary components.
2684
P. Kokkonis, V. Leute / Solid State Ionics 176 (2005) 2681 – 2688
Although the relations between the interdiffusion coefficient for a sublattice and the corresponding diagonal Onsager coefficient cannot exactly be deduced from theoretical reasons, they yield plausible procedures that will allow us to estimate composition dependent diagonal Onsager coefficients needed for the calculation of reaction paths in good accordance with the experimental results.
E ðk; lÞ M12
2.5. The Hessian matrix M To calculate the components of the Hessian matrix we need the mean molar Gibbs energy g (k, l). Analogously to the splitting of g also the M matrix can be split into a standard term, an ideal mixing term and an excess term: g ¼ g 0 þ g m þ g E Y M ¼ M 0 þ Mm þ ME : ð23Þ The standard term M0 can be omitted because all second derivatives of g 0 vanish. The ideal mixing term g m is restricted to the ideal mixing entropy as the ideal mixing enthalpy is zero. h i gm ¼ RT Iln k k ð1 k Þð1k Þ l l ð1 l Þð1lÞ : ð24Þ Thus m ¼ RT I M11
1 1 m m m ; M22 ¼ RT I ; M12 ¼ M21 ¼ 0: k k2 l l2 ð25Þ
The term for the excess energy of a quasiternary a-class system (Mk N1-k )(Xl Y1-l ) was derived in [15]: gE ðk; l Þ ¼ k ð1 k Þ½aðMNÞ þ kIbðMNÞ þ lIcðMNÞ þ kIlIdðMNÞ þ l ð1 l Þ½aðXYÞ þ lIbðXYÞ þ kIcðXYÞ þ lIkIdðXYÞ þ k ð1 l ÞIDR G0 :
ð26Þ
As also shown in [15], the quasiternary interaction parameters a(i), h(i), g(i) and y(i) in Eq. (26) can be calculated from the interaction parameters determined for the four quasibinary subsystems. Treating these binary interaction parameters as depending on temperature means that an excess entropy is considered. The ternary parameter D R G 0 =D R H 0
(a)
T I D R S 0 is the molar standard Gibbs energy for the displacement reaction (1). The components of the excess term ME of the Hessian matrix describing the deviation from ideal behaviour read as follows: E M11 ðk; l Þ ¼ ð2 6k ÞðbðMNÞ þ lIdðMNÞÞ 2ðaðMNÞ þ lIcðMNÞÞ E ðk; l Þ ¼ ð1 2k ÞIcðMNÞ þ 2k 3k 2 IdðMNÞ ¼ M21 þ ð1 2l ÞIcðXYÞ þ 2l 3l2 IdðXYÞ DR G0 E M22 ðk; lÞ ¼ ð2 6l ÞðbðXYÞ þ kIdðXYÞÞ 2ðaðXYÞ þ kIcðXYÞÞ
ð27Þ The ME matrix is symmetrical because of M12E = M21E. The non-diagonal term, being important for the ternary effects, depends on D R G 0 and its dependence on composition is only determined by the parameters c and d. Therefore, in symmetrical systems, where neither the interaction between M and N depends on l (c(MN) = 0), d(MN) = 0) nor the interaction between X and Y on k (c(XY) = 0), d(XY) = 0), the whole non-diagonal term will be reduced to M12 = M21 = D R G 0 [15]. That means that in these symmetrical cases D R G 0 is the only parameter that can cause a thermodynamical coupling between the components of the two sublattices. 3. Simulations of interdiffusion in strict regular quasiternary systems In strict regular systems there is per definition no temperature or composition dependence of the interaction parameters. As shown in Ref. [16] the molar standard Gibbs energy of the reaction MX + NY Y NX + MY (1) is identical with the interaction parameter a(MY), i.e. D R G 0 = a(MY). Thus, the excess part of the M matrix can here be written as: 2aðMNÞ aðMNÞ E : ð28Þ M ¼ aðMYÞ 2aðXYÞ The influence of the thermodynamic properties on the diffusion behaviour of quasiternary systems can best be
(b)
Fig. 1. Reaction paths with a(MY) > 0 for D 2 = 10 I D 1 (a) and D 1 = 10 I D 2 (b).
P. Kokkonis, V. Leute / Solid State Ionics 176 (2005) 2681 – 2688 Table 1 Interaction parameters Quasibinary system i
(Hgk Cd1k )Te
(Hgk Cd1k )Se
Hg(Tel Se1l )
Cd(Tel Se1l )
a(i) / (J/mol)
5900
5900
6300
6300
D R(1)G 0 = a(MY) = 20 500
J/mol;
D R(1)S 0 å 0.
realized by the so-called reaction path C(k(z), l(z)) that can be constructed from the mole fraction profiles k(z) and l(z) [9]. Generally, to explain the diffusion behaviour one uses the projection of the reaction path onto an iso-thermic section l(k) of the phase diagram of the quasiternary system (cf. Fig. 1). The most important feature that determines the shape of a reaction path is a spinodal miscibility gap, because each reaction path touching the spinodal curve crosses the region of instability along a tie line. For a strict regular quasiternary system the isothermal spinodal curve is given by the equation xy 2aðXYÞIx þ 2aðMNÞIy þ 4aðMNÞIaðXYÞ aðMYÞ2 ¼ 0 ð29Þ with x := RT / [k(1 k)] and y := RT / [l(1 l)]. In Ref. [16] it was shown that this means that the shape of the spinodal curve at a given temperature only depends on the magnitude of a(MY), independent of whether the displacement reaction MX + NY Y NX + MY is exo- or endothermic. The shape of the equilibrium curve, however, depends on the sign of a(MY). To show that thermodynamics influences diffusion even in this most simple case, we make very simple assumptions for the L matrix. The kinetic coupling effects shall be negligible which means that the L matrix must be diagonal, i.e. k = 0. Furthermore, the L matrix is taken as independent of composition (k, l) which means that all four diffusion coefficients in Eq. (21) or Eq. (22) must have identical values, i.e. there are only two diffusion coefficients one for the cationic (D 1) and the other for the anionic components (D 2). Fig. 1a and b shows the reaction paths for an endothermic reaction MX + NY Y NX + MY and its reverse exothermic reaction NX + MY Y MX + NY on condition that D 2 = 10 I D 1
2685
(a) or D 1 = 10 I D 2 (b). These reaction paths were simulated by computer calculations using the values for the quasibinary interaction parameters as shown in Table 1. In [4] we have shown that by changing the sign of a(MY) the reaction path for a strict regular system as well as the tie lines and the spinodal curve are reflected at the axis k = 0.5. Therefore only the case a(MY) > 0 is treated here explicitly. Within the spinodal curve the reaction path runs always along the tie lines. Therefore, independent of the diffusion coefficients, the reaction path for an endothermic reaction runs always nearly along the diagonal of the phase square. An exothermic reaction, however, shows quite another behaviour. Here, the reaction path also crosses the spinodal curve via a tie line, but it starts more or less along quasibinary edges until a tie line is reached. If the diffusion coefficient for the chalcogens X and Y is higher than that for the constituents of the metal sublattice M and N (Fig. 1a) the reaction path runs predominantly along the chalcogen edges M(X,Y) and N(X,Y), in the reverse case (Fig. 1b) along the edges (M,N)X and (M,N)Y. Whereas the jump across the miscibility gap is forced by thermodynamics, the shape of the reaction path within the homogeneous region is determined predominantly by kinetics, i.e. by the ratio of the interdiffusion coefficients for the two sublattices. This kinetic effect becomes also apparent near the pure binary components of the endothermic reaction. In [17] we had already demonstrated that a Z-shaped reaction path, as shown in Fig. 1b, is retained even at such high temperatures, where the spinodal miscibility gap does no longer exist. At temperatures above the critical temperature the thermodynamic interaction causes a region of uphill diffusion caused by the off-diagonal part of the Hessian matrix. This behaviour is quite analogous to the influence of the thermodynamic factor on interdiffusion at temperatures above spinodal miscibility gaps in quasibinary systems [3]. Till now the initial conditions (t = 0) for the simulation of the displacement reactions have been the following: ðkV ¼ 0; lV ¼ 1Þ; ðkþV ¼ 1; lþV ¼ 0Þ or ðkV ¼ 1; lV ¼ 0Þ; ðkþV ¼ 0; lþV ¼ 1Þ
Fig. 2. Reaction paths with D 1 = 10 I D2 for a(MY) > 0.
2686
P. Kokkonis, V. Leute / Solid State Ionics 176 (2005) 2681 – 2688
(a)
(c)
(b)
(d)
Fig. 3. Diffusion profiles k(z) and l(z) corresponding to reaction paths selected from Fig. 2 (a(MY) > 0).
But the Hessian matrix M is not influenced by the initial conditions, thus it remains unchanged also for so-called quasibinary initial conditions: ðkV ¼ c; lV ¼ 1Þ; ðkþV ¼ c; lþV ¼ 0Þ; ðFig: 2aÞ or ðkV ¼ 0; lV ¼ cÞ; ðkþV ¼ 1; lþV ¼ cÞ; ðFig: 2bÞ
with
0 V c V 1: The simulations of reaction paths with quasibinary initial conditions, as shown in Fig. 2a and b, yield pronounced deviations from quasibinary sections. These deviations are caused by the ternary interaction parameter a(MY). As long as the reaction path does not intersect the spinodal miscibility gap one obtains reaction paths corresponding to diffusion profiles with regions of uphill diffusion for those components with identical boundary conditions on both ends of the profile; e.g. for interdiffusion of Hg and Cd in Fig. 3a or of Te and Se in Fig. 3c. If the reaction path crosses the spinodal miscibility gap, the diffusion profiles show reverse concentration steps at the phase boundary as for example Hg in Fig. 3b. The comparison of Fig. 3a and b shows that, if the phase diagram were not known, one could hardly distinguish between continuous profiles with uphill diffusion and discontinuous profiles with inverse concentration steps at phase boundaries. Even more complicated are situations where the reaction path cannot cross the spinodal miscibility gap along a tie line because it reaches the spinodal curve near one of the critical points for spinodal demixing as in Fig. 2b for l = 0.75. In such cases one obtains an oscillating behaviour until the reaction path can leave the unstable region (Fig. 3d). 4. The system (Hgk Cd1-k )(Tel Se1-l ) 4.1. Results of computer simulations To show that the kinetic influence on solid state reactions can be reproduced correctly by the computer simulations, diffusion profiles for the displacement reaction HgTe + CdSe Y CdTe + HgSe (1) were simulated with different sets of
Table 2 Experimentally determined interdiffusion coefficients at 723 K [19]; (k, l) are the compositions at which the interdiffusion coefficients were measured, the subscripts Fme_ and Fch_ mean metal or chalcogen saturated ˜ meexp / (cm2/s) ˜ chexp / (cm2/s) Quasibinary subsystem (k, l) D D (Hgk Cd1k )Te (Hgk Cd1k )Te (Hgk Cd1k )Se (Hgk Cd1k )Se Hg(Tel Se1l ) Hg(Tel Se1l ) Cd(Tel Se1l ) Cd(Tel Se1l )
(0.1, (0.9, (0.9, (0.1, (1.0, (1.0, (0.0, (0.0,
1.0) 1.0) 0.0) 0.0) 0.9) 0.1) 0.1) 0.9)
2 I 10 14 3 I 10 12 2 I 10 10 6 I 10 13 3 I 10 14 1.5 I 10 14 2 I 10 17 5 I 10 16
1.4 I 10 13 3 I 10 12 2 I 10 11 3 I 10 13 4 I 10 10 2 I 10 10 – 1 I10 13
diffusion coefficients. The same reaction was already experimentally investigated in dependence on intrinsic doping in [18,19]. The interaction parameters for this system as derived in Ref. [18] are listed in Table 1. We tried to reproduce the diffusion profiles and reaction paths for metal saturated and chalcogen saturated samples at 723 K and for samples with the composition corresponding to minimum vapor pressure at 673 K. To find suitable sets of tracer diffusion coefficients, as defined by Eq. (21) or (22), we used as starting values the interdiffusion coefficients for compositions near the corners of the phase diagram, which were experimentally determined in Ref. [19] for reactions between metal saturated as well as between chalcogen saturated samples (Table 2). The tracer diffusion coefficients estimated by adjusting the calculated diffusion profiles and reaction paths to the experimentally determined data are listed in Table 3. The calculated reaction paths and diffusion profiles are plotted in Figs. 4a –c and 5a –c; they correspond rather good with the experimentally determined ones as shown in [18,19]. It is astonishing that the rather simple assumptions, used for modeling the Onsager matrix, allow to reproduce the experimentally determined diffusion profiles and reaction paths in all essential details. For the simple model used here, all kinetic and thermodynamic data that are needed for the simulation of the diffusion profiles and reaction paths can be taken from experiments on the quasibinary subsystems: & The molar standard Gibbs energy for the displacement reaction D R G 0 follows principally from the standard chemical potentials of the pure binary components. & The ternary interaction parameters (cf. Eq. (26)) can be calculated from the interaction parameters of the quasibinary Table 3 Tracer diffusion coefficients of minor components in the binary compounds determined by simulation; the subscripts Fme_, Fch_ and Fmin_ mean metal saturated, chalcogen saturated or minimal vapor pressure ˜ me / (cm2/s) ˜ ch / (cm2/s) ˜ min / (cm2/s) Tracer diffusion coeficient D D D Hg in CdTe Cd in HgTe Cd in HgSe Hg in CdSe Se in HgTe Te in HgSe Te in CdSe Se in CdTe
1 I10 13 1 I10 12 4 I 10 11 1 I10 12 1 I10 13 1.5 I 10 14 2 I 10 17 1 I10 15
5 I 10 12 3 I 10 12 1 I10 12 3 I 10 13 4 I 10 10 2 I 10 10 1 I10 12 1 I10 12
1 I10 12 1 I10 13 1 I10 13 1 I10 13 5 I 10 12 1 I10 12 1 I10 12 1 I10 12
P. Kokkonis, V. Leute / Solid State Ionics 176 (2005) 2681 – 2688
(a)
2687
(b)
(c)
Fig. 4. Reaction paths for (a) chalcogen and (b) metal saturation and for (c) minimum vapor pressure. The dashed area describes the spinodal miscibility gap with tie lines and the inner curve is the spinodal curve.
subsystems even for nonregular systems. The spinodal miscibility gap– being in accordance with the experimental phase diagram – was calculated with the same set of interaction parameters. & The diffusion coefficients needed for constructing the Onsager matrix according to Eq. (21) or (22) correspond to the extrapolated interdiffusion coefficients of the quasibinary subsystems (Table 2). According to Darken’s equation one obtains by extrapolation of the pure binary components the tracer diffusion coefficients of the minor component, e.g. for the quasibinary subsystem (M,N)X the diffusion coefficients for N in MX and M in NX.
& The coefficient k (Eq. (19)), responsible for kinetically caused ternary diffusion effects, can in most systems with separated sublattices be taken as zero. Corresponding to earlier reported experimental findings [19] the computer simulations yield continuous reaction paths for chalcogen (Fig. 4a) as well as for metal saturated (Fig. 4b) reaction couples. For a reaction couple with the composition of minimal vapor pressure the simulations yield – also corresponding to the experiments [18] – a reaction path crossing the spinodal miscibility gap. 4.2. Uphill diffusion
(a)
(b)
(c)
The diffusion profiles belonging to the continuous reaction paths show distinct zones of uphill diffusion. Starting at the compound HgSe the reaction path runs at first in that direction, where the interdiffusion coefficients are the highest, i.e. for the chalcogen saturated case along the edge Hg(Se,Te) (Fig. 4a) or for the metal saturated case along the edge (Hg, Cd)Se (Fig. 4b). In the metal saturated case the decrease of the Hg concentration is stopped, when the chalcogen exchange starts. In the same measure as the chalcogen exchange increases, the uphill diffusion of Hg gets promoted (cf. Fig. 5b). The essential term for this uphill diffusion is the partial diffusion coefficient D MX of Eq. (12) for M = Hg and X = Te: DHgTe
Fig. 5. Diffusion profiles k(z) and l(z) for (a) chalcogen (t = 1 h) and (b) metal saturation (t = 1 h) and for (c) minimum vapor pressure (t = 5 h).
1 B2 g B2 g LHgHg ¼ þ LHgTe C BxHg BxTe BxT e BxTe ¼
1 LHgHg I aHgTe : C
ð30Þ
2688
P. Kokkonis, V. Leute / Solid State Ionics 176 (2005) 2681 – 2688
For the given system this last equation becomes so simple as L HgTe = 0 because of k = 0 and as in strictly regular systems the off-diagonal terms of the matrix (Eq. (28)) consist only of the ternary interaction parameter a HgTe. Obviously, ternary effects like uphill diffusion can occur without kinetic coupling effects. In most cases uphill diffusion can be explained as caused only by thermodynamics, here for example by a HgTe. The partial diffusion coefficient D HgTe in Eq. (30) has a negative sign because of the positive sign of a HgTe (Table 1). Therefore the Hg flux according to Eq. (31) and Fig. 5b can remain negative even if the gradient flk / flz for this flux becomes negative, provided the slope of the Te profile flk / flz is high enough and also negative. Bk Bl Bk þ DHgTe : ð31Þ JHg ¼ C DHgHg ¼ C D˜ M Bz Bz Bz The sign of the ternary interdiffusion coefficient for the metal sublattice [9], defined by D˜ M : ¼ DHgHg þ DHgTe Bl=Bk;
ð32Þ
depends on the slope of the reaction path, fll / flk. An uphill diffusion of the cationic components demands (Eq. (31)) a ˜ M and thus, because of D HgHg > 0 (Eq. negative value for D (18)), it demands opposite signs for DHgTe and fll / flk. Furthermore, the magnitude of the second term in Eq. (32) must be high enough to over-compensate the first term. This can be easier attained the higher the slope of the reaction path. The discussion of uphill diffusion for the chalcogen saturated case has to proceed analogously. 5. Conclusions The calculations of interdiffusion profiles and reaction paths for solid state reactions in a-class systems have demonstrated that the experimentally observed phenomena can be explained by a rather simple model. Whether the reaction path has to cross miscibility gaps or runs along continuous solid solutions depends on the Hessian matrix. Especially ternary effects like uphill diffusion can be understood as caused predominantly by thermodynamic coupling effects initiated by the interaction parameter a(MY) being identical with D R G 0 of the displacement reaction MX + NY Y NX + MY. The influences of intrinsic doping on diffusion are described by the Onsager matrix. In those calculations where the reaction path crosses the spinodal miscibility gap a problem arises, because the region between the spinodal curve and the equilibrium curve of the spinodal miscibility gap is still stable with regard to diffusion [20]. Therefore, a reaction path that is calculated with the described simulation model does not become unstable before it reaches the spinodal curve. Under conditions of local equilibrium, however, the reaction path should end already at the equilibrium curve. To solve this problem the precipitation processes within the metastable region would have to be considered, i.e. nucleation kinetics would have to be
included into the diffusion algorithms, a nearly unsolvable task. But compared with experiments, this problem is not so serious at a glance. A diffusion profile crossing a phase boundary, determined experimentally by an electron microprobe, can hardly be distinguished from a calculated diffusion profile, because of the sharp concentration changes near the interface. Calculated profiles yield in this small region only a few calculated points and measurements with a microprobe do also not yield sharp concentration steps, but a more or less continuous course, because of the restricted lateral resolution of the measured X-ray signals. Thus, considering that one cannot even be sure that during diffusion processes local equilibrium at the interface is really established, the calculated profiles correspond unexpectedly well with the experimental results. All things considered, the simulations of reaction paths and concentration profiles, on the basis of the described model, allow us to understand ternary effects like uphill diffusion and to estimate sets of diffusion coefficients by which the essentials of measured concentration profiles can be reproduced. (The computer programs for the simulation processes can be downloaded from the Internet [21].) References [1] P. Glansdorff, I. Prigogine, Thermodynamic Theory of Structure, Stability and Fluctuations, Wiley-Interscience, London, 1971, p. 30. [2] A.R. Allnatt, A.B. Lidiard, Atomic Transport in Solids, Cambridge Univ. Press, Cambridge, 1993. [3] V. Leute, Solid State Ion. 164 (2003) 159. [4] A. Polyzois, Kokkonis, Die Rolle der Thermodynamik fu¨r die Deutung von Kopplungseffekten bei mehrkomponentigen Diffusionsprozessen Thesis (1997) Mu¨nster. [5] J.S. Kirkaldy, D.J. Young, Diffusion in the Condensed State, The Institute of Metals, London, p. 154f. ˚ gren, J. Phys. Chem. Solids 43 (1982) 421. [6] J. A [7] W. Schottky, in: W. Schottky (Ed.), Halbleiterprobleme, vol. 4, Vieweg, Braunschweig, 1958, p. 235. [8] F.A. Kro¨ger, The Chemistry of Imperfect Crystals, Vol. 1-3, North Holland, Amsterdam, 1974. [9] V. Leute, Solid State Ion. 17 (1985) 185. [10] A.R. Allnatt, A.B. Lidiard, Atomic Transport in Solids, Cambridge Univ. Press, Cambridge, 1993, p. 188. [11] S.R. de Groot, P. Mazur, Grundlagen der Thermodynamik Irreversibler Prozesse, Bibliograph. Inst., Mannheim, 1969, p. 35. [12] L. Onsager, Phys. Rev. 37 (1931) 405. [13] A.R. Allnatt, A.B. Lidiard, Atomic Transport in Solids, Cambridge Univ. Press, Cambridge, 1993, p. 197. [14] J.S. Kirkaldy, D.J. Young, Diffusion in the Condensed State, The Inst. of Metals, London, 1987, p. 142. [15] V. Leute, Ber. Bunsenges. Phys. Chem. 93 (1989) 7. [16] V. Leute, Ber. Bunsenges. Phys. Chem. 80 (1976) 866. [17] P.A. Kokkonis, V. Leute, Defect Diffus. Forum 143/147 (1997) 1159. [18] V. Leute, H.M. Schmidtke, W. Stratmann, Ber. Bunsenges. Phys. Chem. 86 (1982) 732. [19] V. Leute, H.M. Schmidtke, W. Stratmann, A. Kalb, Ber. Bunsenges. Phys. Chem. 87 (1983) 483. [20] I. Prigogine, R. Defay, Chemische Thermodynamik, VEB Deutscher Verlag fu¨r Grundstoffindustrie, Leipzig, 1962, p. 243. [21] http://www.uni-muenster.de/Chemie.pc/Leute/Welcome.html.