5516
J. Phys. Chem. 1996, 100, 5516-5524
Green-Kubo Formulas for Mutual Diffusion Coefficients in Multicomponent Systems Yanhua Zhou* and Gregory H. Miller Department of Geophysical Sciences, UniVersity of Chicago, 5734 South Ellis AVenue, Chicago, Illinois 60637 ReceiVed: NoVember 16, 1995X
Green-Kubo relations for the mutual diffusion coefficients of an n-component system are derived from hydrodynamic equations, linear response theory, and the statistical theory of fluctuations. The mutual diffusion coefficients, which form an (n - 1)-dimensional matrix, are formulated as a product of a kinetic matrix and a thermodynamic matrix. The kinetic matrix consists of time integrals of diffusion flux correlation functions, i.e., the Onsager phenomenological coefficients. The thermodynamic matrix is determined by spatial integrals of radial distribution functions. The formulas are derived initially in the barycentric reference frame and then transformed to the volume-fixed reference frame, which is more appropriate to laboratory studies. Explicit expressions for binary and ternary mixtures are presented, from which various empirical relations between the mutual diffusion coefficients and self-diffusion coefficients are obtained using particular assumptions regarding the velocity cross-correlation functions. Assuming zero values for certain combinations of the cross-correlation functions yields the Darken or Hartley-Crank equation for binary systems and Cooper’s diffusion model for ternary systems, while assigning negative values to these combinations results in Lasaga’s relation, which is a special case of the so-called common force diffusion model.
1. Introduction The statistical-mechanical theory of irreversible processes establishes a connection between microscopic time-correlation functions and macroscopic transport coefficients, regardless of whether the deviation from equilibrium is caused by externally imposed perturbations or naturally-occurring fluctuations.1-4 The mathematical formulation of this fluctuation-dissipation theorem is an equation for the phenomenological coefficient of a transport process in terms of the Fourier component of timedependent fluctuations of the appropriate microscopic dynamical variables, known as Green-Kubo formulas.4 The essence of the theorem is that one may determine transport properties of a system from the fluctuations around the equilibrium state: macroscopic transport need not occur.5 It provides a sound basis for deducing transport coefficients from equilibrium molecular dynamics simulations, which have greatly enhanced our understanding of dynamical processes in condensed matter (e.g., ref 6). Diffusion as a mass transport mechanism has been the subject of numerous theoretical, experimental, and computer simulation studies. Its complexity goes far beyond the apparent simplicity of Fick’s law. It is particularly complicated in multicomponent systems, where unusual phenomena such as uphill diffusion7 and double diffusive convection8 are observed. To describe the diffusion processes in a multicomponent system three types of diffusion coefficients have been defined: self-diffusion and trace and mutual diffusion coefficients (the last are also called interdiffusion coefficients9 or chemical diffusion coefficients10). Thorough discussions of these quantities and their relationships to laboratory measurements are given in refs 9 and 10. Although diffusion in binary systems has been extensively studied, the study of diffusion in ternary and higher-order systems is still in a comparatively exploratory stage. Mutual diffusion coefficients have been measured in only a few multicomponent systems,11-13 and to our knowledge no computer simulation studies have calculated mutual diffusion coefficients in any system of more than two components. The X
Abstract published in AdVance ACS Abstracts, March 1, 1996.
0022-3654/96/20100-5516$12.00/0
related statistical theory is also incomplete: whereas expressions for the binary mutual diffusion coefficients are well-known,3,14,15 the analogous expressions for the mutual diffusion coefficients in ternary and higher-order systems have not been published. One objective of this study is to derive the Green-Kubo formulas for mutual diffusivities of generic multicomponent systems. This study is principally motivated by an attempt to understand commonly-used empirical diffusion models from a molecular point of view. These empirical diffusion models, which appear in the chemical, metallurgical, and geochemical literature,16-20 give relations for the mutual diffusion coefficients in terms of self-diffusion coefficients and thermodynamical parameters such as partial molar volumes and mole fractions. These models are of interest because they bridge the gap between simple end member laboratory experiments and more complicated natural systems.21,22 For example, in natural geological systems dozens of major, minor, and trace elements interdiffusion during processes of melting, crystallization, melt percolation, magma mixing, and fluid infiltration. It is unfeasible to experimentally determine all (n - 1)2 mutual diffusion coefficients for every system of interest. Given the empirical diffusion models, however, one may estimate the mutual diffusion coefficients of the natural system from a knowledge of the appropriate self-diffusion coefficients. In this study we will show how the empirical diffusion models may be deduced as special cases of the complete Green-Kubo mutual diffusion coefficient formulas. In the process, the importance of velocity cross-correlation functions will be demonstrated. We begin the derivation with the hydrodynamic and linear phenomenological equations. For ease and clarity of presentation we consider diffusion in an isothermal and isobaric system, driven only by chemical-potential gradients. In the linearized theory the mutual diffusion coefficient expressions we obtain are unaffected by this simplification. A more general treatment would obtain coefficients for the Soret and Dufour effects, as well as other phenomenological coefficients. In section 2 the related fundamental equations are first outlined. Then, the Green-Kubo formulas for mutual diffusion coefficients in an © 1996 American Chemical Society
Green-Kubo Formulas for Diffusion Coefficients
J. Phys. Chem., Vol. 100, No. 13, 1996 5517
n-component system are derived in the barycentric reference frame. A transformation of the results to the volume-fixed reference frame is presented at the end of the section. In section 3 the correlation functions involved in the formulas are expressed in terms of appropriate microscopic particle variables and functions. Complete microscopic expressions for the mutual diffusion coefficients of binary and ternary systems are presented in section 4. Section 5 shows how the relations between mutual diffusion and self-diffusion coefficients proposed in empirical diffusion models may be obtained using specific simplifying assumptions regarding the velocity cross-correlation functions.
where
GRβ )
Fm
∂cR(b, r t) ∂t
r t) ) 0 + ∇‚J BR(b,
r t) ) B J R(b,
1
r t) ∂cR(b, ∂t
∑LRβ∇µβm(b,r t)
(2.2)
where T is the equilibrium temperature, LRβ are the Onsager r t) is given by µmR ) phenomenological coefficients, and µmR (b, µR/mR - µn/mn, with mR and µR the molar mass and molar chemical potential of component R. Since mutual diffusion coefficients are defined in terms of concentration gradients, rather than chemical-potential gradients, we may alternatively write the linear response equation (2.2) in the following form:15
∂t
c˜ Rk (ω) ) ∫0 dt ∫V db r exp(-ik B‚b r + iωt) cR(b, r t) ∞
where DRβ are the mutual diffusion coefficients, which form an (n - 1)-dimensional matrix in an n-component system. While eq 2.2 may be considered more fundamental than (2.3), the latter has greater practical value since concentration gradients are more readily measured. The relation between DRβ and LRβ is found immediately by equating (2.2) and (2.3). In matrix notation the relation takes the form
D)LG
(2.4)
(2.7)
(2.8)
where we have used z ) -iω for the Laplace transformation argument and it may have a small negative real part to provide necessary integral convergence. From the spatial Fourier transform of eq 2.1, the continuity equation in k-space is
∂cRk (t) Fm + ik B‚J BRk (t) ) 0 ∂t
(2.9)
The diffusion equation (2.7) in (ω, k)-space is
b ck c˜ k(ω) ) ∆-1 b
(2.10)
∆ ) -iωI + k2D
(2.11)
where
ck are the (n with I the identity matrix. In eq 2.10 b c˜ k(ω) and b - 1)-dimensional column vectors whose Rth elements are the b, t) and the Fourier transform Fourier-Laplace transform of cR(r b, 0), respectively. of cR(r 2.2. Statistical-Mechanical-Correlation Functions. From (2.10) an (n - 1) × (n - 1)-matrix equation composed of concentration-correlation functions can be constructed
∆-1 ) ΩC-1
(2.3)
β)1
c (b, r t) ) D ∇2 b
We define the combined Fourier and Laplace transforms of b, t) to be c˜ Rk (ω), cR(r
n-1
r t) ) -Fm ∑ DRβ ∇cβ(b, r t) B J R(b,
(2.6)
β)1
∂c b(b, r t)
n-1
Tβ)1
n-1
r t) ) ∑ DRβ ∇2 cβ(b,
For convenience we introduce an (n - 1)-dimensional vector b c(r b, t) consisting of the set of variables c1(r b, t), c2(r b, t), ..., cn-1(r b, t). It then follows from (2.6) that
(2.1)
where Fm is the mass density of the system. There are (n - 1) independent equations with the form of (2.1). The relation n B JR(r b, t) ) 0 completes this system of equations and ∑R)1 defines the barycentric reference frame, a reference frame appropriate for molecular dynamics simulations and the one we use now. Without sacrificing generality we choose the nth component as the dependent variable. To complete the hydrodynamic description the Navier-Stokes equations and the equation of energy conservation should also be included, but since we are only interested in the isothermal and isobaric systems, they are ignored here. In the linearized theory irreversible flows are linear functions of “thermodynamic forces”. Having assumed that the local temperature and pressure are maintained at equilibrium, the only thermodynamic forces remaining in the system are spatial gradients of chemical potentials. We may write the linear phenomenological law for the diffusion fluxes as23
(2.5) T,P,cγ*β
The matrix L will be shown in the following sections to depend on time integrals of diffusion flux correlation functions, and the strictly thermodynamic matrix G will be related to spatial integrals of radial distribution functions. A form of Fick’s second law is obtained by combining eqs 2.1 and 2.3:
2. General Expressions for Mutual Diffusion Coefficients 2.1. Hydrodynamic and Phenomenological Equations. b, t) denote the mass fraction of component R at Let cR(r coordinate b r and time t in an n-component system. The rate of b, t), by concentration change is linked to the diffusion flux, B JR(r the conservation of mass,23 i.e.,
( )
∂µmR Fm T ∂cβ 1
(2.12)
where the matrices Ω and C consist of time-dependent and static concentration-correlation functions, respectively, given by β 〉 ΩRβ ) 〈c˜ Rk (ω)c-k
(2.13)
β 〉 CRβ ) 〈cRk c-k
(2.14)
and
The symbol 〈 〉 represents an equilibrium ensemble average over b, 0). the fluctuations in the initial values cR(r
5518 J. Phys. Chem., Vol. 100, No. 13, 1996
Zhou and Miller
In general we can write the inverse of a matrix M explicitly as [M-1]βR ) |M|Rβ/|M|, where |M|Rβ is the cofactor of MRβ in the determinant |M|. In this notion eq 2.12 can also be written
|∆|βR |∆|
n-1
1
)
∑ ΩRγ|C|βγ
(2.15)
|C| γ)1
Now taking the inverse Laplace transform of eq 2.15 gives
|∆|βR - δRβ ) (-iω) |∆| 1 ∞d
∫ |C| 0
dt
D11 )
1 2 3Fm 〈|c1|2〉
γ 〉|C|βγ] exp(iωt) dt (2.16) [ ∑ 〈cRk (t)c-k γ)1
[
|∆|βR
lim (-iω)
|∆|
]
- δRβ )
DRβk2 (2.17)
iω
and
[
D11 D21
D12 D22
]
ωf0 kf0
iω
d
∫0
∞
k2|C|
1 2 3Fm|C0|
[
∫0∞〈JB1(t)‚JB1〉 dt ∫0∞〈JB1(t)‚JB2〉 dt ∫0∞〈JB2(t)‚JB1〉 dt ∫0∞〈JB2(t)‚JB2〉 dt
[
dt
n-1
γ 〉|C|βγ] exp(iωt) dt [ ∑ 〈cRk (t)c-k γ)1
(2.18)
In eq 2.18 it is essential to take the limit in the order k f 0 before ω f 0. The time-dependent concentration-correlation functions in eq 2.18 can be replaced by the corresponding diffusion fluxcorrelation functions. This is accomplished by using eq 2.9 and the time-translational invariance property of correlation 1 1 〉 ) 〈J B1k c-k ( - t)〉. An integration by functions, i.e., 〈J B1k (t) c-k parts of eq 2.18 gives the result
DRβ ) lim lim ωf0 kf0
1
∫0
∞
2 3Fm |C|
n-1
R [ ∑ 〈J BRk (t)‚J B-k (t)〉|C|βγ] × γ)1
exp(iωt) dt )
2 3Fm |C0|
∫0∞ [∑〈JBR(t)‚JBγ〉|C0|βγ] dt
where the matrix C0 is given by [C0]Rβ ) 〈cR cβ〉 and is obviously symmetric. The subscript k ) 0 for B J R0 and cR0 has been dropped to simplify the notation, and the superscripted component indices in the expression remind us that these quantities are Fourier transformed. Putting eq 2.19 in the same form as eq 2.4, we have
D)LG where L and G are given by
LRβ )
1 ∫∞ 〈JBR(t)‚JBβ〉 dt 3kBV 0
GRβ )
|C0|Rβ
2 Fm |C0|
TRβ ) Fm
[
(2.24)
(
hβ V hn V δRβ + FR mR mn mβ
()
DV ) TD
respectively, with V the total volume of the system and kB the
(2.23)
)]
(2.25)
with V h R and FR the partial molar volume and molar density of component R. In addition, if we change mass fraction gradients in eq 2.3 to molar density gradients, which in some cases are more appropriate to laboratory experiments, and substitute in eqs 2.4 and 2.21, we may further write eq 2.24 as
(2.20)
(2.21)
]
where T is the (n - 1)-dimensional reference frame transformation matrix, which in this case is given by
and
k BV
×
〈|c2|2〉 -〈c1 c2〉 -〈c1 c2〉 〈|c1|2〉
DV ) T D
(2.19)
γ)1
]
where |C0| ) 〈|c1|2〉〈|c2|2〉 - (〈c1 c2〉)2. The binary result of eq 2.22 is well-known,3,14,15 but the ternary result of eq 2.23 has not appeared previously. From the ternary formulas we can B1(t)‚J B1〉 ) see that D12 is generally not equal to D21 unless 〈J 2 2 1 2 2 2 B 〉 and 〈|c | 〉 ) 〈|c | 〉. 〈J B (t)‚J 2.3. Reference Frame Transformation. The formulas derived above are valid in the barycentric reference frame. We adopt this reference frame for its natural suitability to the microscopic description of systems, and molecular dynamics simulations in particular. However, in most laboratory diffusion experiments the volume-fixed reference frame is more appropriate.25,26 In the remainder of this section we transform the above results to this reference frame. As shown by De Groot and Mazur23 the diffusivity matrix in the volume-fixed reference frame, DV, is related to that in the barycentric reference frame, D, by the relation
n-1
1
(2.22)
)
from which the following expression for DRβ is obtained immediately,
DRβ ) lim lim
∫0∞ 〈JB1(t)‚JB1〉 dt
n-1
In the zero wavenumber limit the left-hand side of the above equation is
kf0
Boltzmann constant. In eq 2.20 the matrix L is expressed in terms of time integrals of diffusion flux-correlation functions. It is symmetric when 〈J BR(t)‚J Bβ〉 ) 〈J Bβ(t)‚J BR〉, which follows from the principle of microscopic reversibility.24 As will be shown in the following section, expression 2.21 for the matrix G is identical to that given by (2.5). From eqs 2.19-2.21 we can easily obtain the results for binary (n ) 2) and ternary (n ) 3) systems:
)
kB V 2 Fm
∂c b
∂F B
T L C-1 0
() ∂c b
∂F B
(2.26)
The molar density vector B F is of (n - 1)-dimension and its nth
Green-Kubo Formulas for Diffusion Coefficients
J. Phys. Chem., Vol. 100, No. 13, 1996 5519
n V h RFR ) 1. component can be obtained from the relation ∑R)1 The last term in eq 2.26 is given by
( ) ∂cR ∂Fβ
) T,P,Fγ*β
[
(
)]
V hβ 1 m δ + cR mn - mβ Fm R Rβ V hn
3. Related Microscopic Expressions 3.1. Diffusion Flux B Jr(t). In the expression (2.20) for the k ) 0 component of the Fourier transform matrix L, B JR(t) is the B b, t). A microscopic expression for this of the diffusion flux B JR(r term can be found by considering concentration fluctuations in the system. In an n-component system with fixed total volume, mass fraction fluctuations may be written in terms of number density fluctuations,
mR
B J R(t) ) mRBj R(t)
(2.27)
From (2.25-2.27) we find that the relationship for the mutual diffusivities between the two reference frames is simple for binary systems, i.e., D11 ) DV11, but relatively complicated for higher-order systems.
δcR )
n In the barycentric reference frame where ∑R)1 mRBjR(t) ) 0 the diffusion fluxes take a particularly simple form,
However, in the volume-fixed reference frame where instead n V h RBjR(t) ) 0, there is no such simplification, and the ∑R)1 general result eq 3.5 must be used. 3.2. Matrix C0. The microscopic expression for 〈cR cβ〉 can be derived as follows. From (3.1) we know that 〈cR cβ〉 are related to 〈FR Fβ〉, which, in turn, can be written in terms of 〈NRNβ〉 using eq 3.4a. Furthermore, following Kirkwood and Buff,27 〈NRNβ〉 are linked to radial distribution functions gRβ(r) in the canonical ensemble:
N)VB
∑ (mj δRβ - xRmβ)δFβ
(3.1)
where m j ) Fm/F is the average mass per particle and xR ) FR/F the mole fraction of component R. From the microscopic perspective we can define number densities and particle currents in terms of individual particles’ coordinates and velocities. The Fourier transform of the number b, t) and BjR(r b, density and particle current of component R, FR(r t), can be written as
(3.7)
where the matrices N and B are given by
n
Fmm j β)1
(3.6)
NRβ ) 〈NRNβ〉 - 〈NR〉〈Nβ〉
(3.8)
BRβ ) FRδRβ + FRFβΓRβ
(3.9)
ΓRβ ) 4π ∫0 r2[gRβ(r) - 1] dr
(3.10)
and
with ∞
NR
FRk (t)
) ∑ exp[-ik B‚b r i(t)]
(3.2a)
i)1
and NR
u i(t) exp[-ik B‚b r i(t)] Bj Rk (t) ) ∑ b
(3.2b)
i)1
where NR is the number of particles of component R, and b ri(t) and b ui(t) are the coordinate and velocity vectors of particle i of component R. In k-space the conservation of component R is15
∂FRk (t) + ik B‚Bj Rk (t) ) 0 ∂t
(3.3)
It can be easily seen that in the limit k f 0 eqs 3.2 reduce to R
F (t) ) NR
(3.4b)
i)1
In combination, the above equations and eq 2.9 give an expression for the diffusion fluxes in terms of the particle currents or the particle velocities,
mR
(3.11)
with
ARβ )
NR
B J R(t) )
A ) (V B)-1
(3.4a)
and
u i(t) Bj R(t) ) ∑ b
Combining these relations we can obtain expressions for C0 in terms of the microscopic structural functions ΓRβ. The final results are given in Appendix B. Note that 〈cR cβ〉 is an (n 1)-dimensional matrix, whereas the matrices 〈FR Fβ〉, N, and B are n-dimensional. Finally, to relate the matrix C0 to the derivatives of chemical potential with respect to mass fraction, we need relations between the chemical potential derivatives and the radial distribution functions. Kirkwood and Buff27 showed that the derivatives of chemical potential with respect to particle number in the canonical ensemble can be related to the radial distribution functions by
( )
∂µR kB T ∂Nβ 1
(3.12) T,V,Nγ*β
From all the relations discussed in this subsection, and through some lengthy but straightforward algebraic manipulation, we can show that
( )
Fm ∂µmR |C0|Rβ ) |C0| V kBT ∂cβ
(3.13)
T,P,cγ*β
n
(m j δRβ - xRmβ)Bj β(t) ∑ m j β)1
(3.5)
Substitution of this result into (2.21) gives eq 2.5 for the thermodynamic matrix G.
5520 J. Phys. Chem., Vol. 100, No. 13, 1996
Zhou and Miller
4. Matrix DV for Binary and Ternary Systems
C), we can write the ternary diffusivity matrix DV as
With the help of (3.13) the diffusivity matrix equation, eq 2.26, may be written as
DV )
1 FmT
TL
( )() ∂µ b
∂x b
3V x3
∂x b
m
∂F B
T,P
(4.1)
() ∂Fβ
) T,P,Fγ*β
1
[ ( )]
F
δRβ + xR
V hβ V hn
-1
where Bjr(t)
1
1
1
2
2
1
2
2
m1m2 r Bj (t) m j
Bj r(t) ) x2Bj 1(t) - x1Bj 2(t)
Substituting this and the result for (∂µ1/∂x1) from Appendix C into eq 4.1, we obtain the binary mutual diffusion coefficient:
Q ∫∞ 〈Bj r(t)‚Bj r〉 dt 3Nx1x2 0
(4.5)
where the thermodynamic factor15 Q may be expressed in terms of the three radial distribution functions (see Appendix C). 4.2. Ternary Systems. The two independent diffusion fluxes for a ternary system are
B J 1(t) )
m1 [m3Bj r1(t) + m2Bj r3(t)] m j
(4.6a)
B J 2(t) )
m2 [m3Bj r2(t) - m1Bj r3(t)] m j
(4.6b)
and
where Bjri(t) are the relative particle currents:
Bj r1(t) ) x3Bj 1(t) - x1Bj 3(t) Bj r2(t) ) x3Bj 2(t) - x2Bj 3(t)
(4.7a) (4.7b)
Q11/x1 Q12/x2 Q21/x1 Q22/x2
]( ) ∂x b
∂F B
(4.8)
T,P
(4.9a)
h 3Bj r2(t) - V h 1Bj r3(t) Bj V2(t) ) V
(4.9b)
and
In eq 4.8 Qij are the thermodynamic factors of the ternary system which may be expressed in terms of radial distribution functions. Of the four Q-factors only three are independent (cf. Appendix C). This is consistent with the fact that there are only three independent chemical-potential derivatives with respect to particle number for a ternary isothermal and isobaric system (because of the Gibbs-Duhem and Maxwell relations). 5. Relations between Mutual Diffusion and Self-Diffusion Coefficients The self-diffusion coefficient of component R in any system is15
(4.3)
(4.4)
×
h 3Bj r1(t) + V h 2Bj r3(t) Bj V1(t) ) V
DR )
is the relative particle current,
D11 ) DV11 )
[
]
where the new relative particle currents BjVi(t) are the relative particle currents Bjri(t) weighted by the partial molar volumes,
(4.2)
In eq 4.1 the diffusivity matrix DV is expressed entirely in terms of microscopic variables. The matrix T and the last matrix are given by eqs 2.25 and 4.2, respectively. From eq 2.20 the matrix L is determined by the time-correlation functions of diffusion fluxes, which are related to particle velocities and mole fractions via eqs 3.5 and 3.4b. The matrix (∂µ b/∂x b)T,P, describing the thermodynamic behavior of the mixture, is a function of integrals of the radial distribution functions ΓRβ, as discussed in the preceding section. Explicit expressions for (∂µ b/∂x b)T,P in binary and ternary systems are presented in Appendix C. These equations comprise a prescription for the computation of multicomponent mutual diffusion coefficients in molecular dynamics simulations. 4.1. Binary Systems. From (3.5) we may write the only independent diffusion flux in binary mixtures as
B J 1(t) )
[
∫0∞〈Bj V (t)‚Bj r 〉 dt ∫02〈Bj V (t)‚Bj r 〉 dt ∫0∞〈Bj V (t)‚Bj r 〉 dt ∫0∞〈Bj V (t)‚Bj r 〉 dt
T,P
b/∂F B)T,P ) (∂x b/ where we have used the relation (∂x b/∂c b)T,P (∂c ∂F B)T,P, which is given by
∂xR
F
NR
1 3NR
∫0 ∑ i)1
∞
〈u bi(t)‚u bi〉 dt
(5.1)
Unlike mutual diffusion coefficients, self-diffusion coefficients are independent of reference frame. Note that two differences exist between eq 5.1 and eq 4.5 or eq 4.8. First, the self-diffusion coefficients are independent of the chemical potential derivatives whereas the mutual diffusion coefficients depend on (∂µ b/∂x b)T,P through the thermodynamic Q-factors. Second, the velocity correlation functions in the mutual diffusion coefficient are collective: the velocities of every particle in the system contribute to the correlation functions, whereas the selfdiffusion coefficient depends on velocity autocorrelation functions of individual particles. These differences explain why in numerical simulations it is difficult to achieve the same precision in mutual diffusion coefficients as in self-diffusion coefficients.28-30 In the following we will show how the mutual diffusion coefficients can be related to the self-diffusion coefficients under different simplifying assumptions. 5.1. Binary Case. If we group separately the velocity autocorrelation and cross-correlation functions that appear in the binary mutual diffusion coefficient and use the f-notation of McCall and Douglass18,31 for the cross-correlation functions, eq 4.5 can be written as
[
(
DV11 ) Q x2D1 + x1D2 + x1x2
f11 x21
+
f22 x22
f12 -2 x1x2
)]
(5.2)
where the f-factors are given by
and
Bj r3(t) ) x2Bj 1(t) - x1Bj 2(t)
(4.7c)
Combining these and expressions for (∂µR/∂xβ)T,P (see Appendix
fRβ )
1
NR Nβ
∑∑∫ 3N i)1 j*i 0
∞
〈u bi(t)‚u bj〉 dt
(5.3)
Green-Kubo Formulas for Diffusion Coefficients
J. Phys. Chem., Vol. 100, No. 13, 1996 5521
If we assume that the velocity cross-correlations satisfy
f11 x21
+
f22 x22
f12 -2 )0 x1x2
(5.4)
a simple relation connecting the mutual diffusion and selfdiffusion coefficients in a binary system follows immediately,
DV11
) Q(x2D1 + x1D2)
(5.5)
Relation 5.5 is often referred to as Darken’s equation16,32 or the Hartley-Crank relation.17,18 For a thermodynamically ideal mixture, the three radial distribution functions are identical, Γ11 ) Γ22 ) Γ12,17 and the thermodynamic factor Q is equal to 1 (cf. eq C1). Binary mixtures of Lennard-Jones fluids are approximately ideal in this sense.15 Equation 5.5 describes diffusion in argon-krypton mixtures of various compositions;28,29,33,34 but a positive deviation from eq 5.5 is found for highly dissymmetric mixtures of Lennard-Jones fluids.35 Here we refer to the case DV11|exp > Q(x2D1 + x1D2) as a positive deviation. For other more complex binary nonelectrolyte solutions, laboratory experimental data reveal both positive and negative deviations from eq 5.5.18 However, (5.5) is approximately satisfied when such solutions exhibit nearly ideal mixing, e.g., a mixture of cyclooctane and cyclopentane.26 As pointed out by Mills et al.26 in a recent study of binary nonelectrolyte mixtures, completely random motions between two particles may result in their velocity cross-correlation functions being approximately zero. Positive and negative velocity crosscorrelation functions may result from attractive and repulsive interparticle interactions, respectively. A similarity exists between eq 5.5 and the Nernst-Einstein relation, an empirical relation that links the electrical conductivity of an ionic system to the self-diffusion coefficients of the cations and anions in the system.36 The electrical conductivity of a molten alkali-metal salt such as NaCl is proportional to the same collective velocity-correlation function as that appearing in the binary mutual diffusion coefficient. The only difference is that the mole fractions, x1 and x2, in the ionic system are no longer variables but are fixed by the constraint of overall charge neutrality. For such an ionic system the crosscorrelation functions, calculated in a molecular dynamics simulation study,36 are such that a negative deviation from eq 5.4 is exhibited. In a study of the diffusion in silicate melts Lasaga20 proposed
DV11 )
QD1D2 x1D1 + x2D2
f12 (D1 - D2) + 2 -2 )2 x x x x1 x2 1 2 1D1 + x2D2 f22
1 × 3Nx3
[
∫0∞〈[Bj r (t) + Bj r (t)]‚Bj r 〉 dt ∫0∞〈[Bj r (t) + Bj r (t)]‚Bj r 〉 dt ∫0∞〈[Bj r (t) - Bj r (t)]‚Bj r 〉 dt ∫0∞〈[Bj r (t) - Bj r (t)]‚Bj r 〉 dt
(5.7)
In this case the net contribution of the cross-correlation functions to the mutual diffusivity is also negative, consistent with the electrical study discussed above. It is easily seen that the mutual diffusion coefficient (5.6) is always less than or equal to that given by (5.5). 5.2. Ternary Case. To simplify the mutual diffusion coefficients in a three-component system let us first assume that the partial molar volumes of all the components in the system
3
1
1
3
2
2
3
1
2
3
2
fRR x2R
[
Q11/x1 Q12/x2 Q21/x1 Q22/x2
]
]
×
(5.8)
+
fββ xβ2
fRβ -2 )0 xRxβ
(5.9a)
and
fRR x2R
+
fβγ fRβ fRγ )0 xβxγ xRxβ xRxγ
(5.9b)
Given these equations, we can write eq 5.8 as
[
]
x1[(x2 + x3) D1 + x1D3] x1x2(D3 - D2) × x1x2(D3 - D1) x2[(x1 + x3) D2 + x2D3] Q11/x1 Q12/x2 Q21/x1 Q22/x2 (5.10)
[
]
which is the analogous relation to (5.5) for ternary systems, and often referred to as Cooper’s diffusion model.19,21 Similarly, the generalized version of assumption 5.6 is
fRR
+
x2R
fββ
fRβ 1 -2 ) - (DR - Dβ)2 2 xRxβ D h xβ
(5.11a)
and
fRR
+
x2R
fβγ xβxγ
-
fRβ
-
xRxβ
fRγ
)
xRxγ
1 - (D2R + DβDγ - DRDβ - DRDγ) (5.11b) D h where D h is the molar average self-diffusion coefficient of the system, D h ) x1D1 + x2D2 + x3D3. From (5.11) Lasaga’s relation for the ternary case is obtained
1 D h
2
1
A generalization of the assumption (5.4) to the ternary case yields two equations for the cross-correlation functions,
(5.6)
where the denominator is the molar average self-diffusion coefficient of the system. Comparing with eq 5.2 we find that eq 5.6 results from the assumption
f11
are approximately equal, i.e., V h1 = V h2 = V h 3 = 1/F. With this assumption the matrix (∂x b/∂F B) reduces to I/F, and the volumeweighted relative particle currents BjVi(t) are also simplified. The diffusivity matrix DV can then be written as
[
×
]
x1D1 [x2D2 + (x1 + x3) D3] x1x2D1 (D3 - D2) × x1x2D2 (D3 - D1) x2D2[x1D1 + (x2 + x3) D3] Q11/x1 Q12/x2 Q21/x1 Q22/x2 (5.12)
[
]
This relation is a version of the “common force model” described in a recent experimental study of diffusion in ternary silicate melts.22 Now if we additionally assume that the ternary mixture is thermodynamically ideal, all six radial distribution functions of
5522 J. Phys. Chem., Vol. 100, No. 13, 1996
Zhou and Miller
the system are identical and the thermodynamic factors Qij are equal to δij (Appendix C). Applying these further simplifications to eqs 5.10 and 5.12, we have
[
]
DV11 DV12 ) DV21 DV22
[
[(x2 + x3) D1 + x1D3] x1(D3 - D2) x2(D3 - D1) [(x1 + x3) D2 + x2D3]
]
(5.13)
for Cooper’s model, and
[ [
DV11 DV21
DV12 DV22
]
)
1 D h
6. Conclusions
×
D1[x2D2 + (x1 + x3) D3] x1D1 (D3 - D2) x2D2 (D3 - D1) D2[x1D1 + (x2 + x3) D3]
]
(5.14)
for the common force model. It can be seen from eqs 5.13 and 5.14 that the model-predicted mutual diffusion coefficients satisfy the stability criterion,37 det(DV) ) DV11DV22 - DV12DV21 > 0. In addition, the predicted diagonal diffusion coefficients DV11 and DV22 are always positive, while the off-diagonal diffusion coefficients DV12 and DV21 may be of any sign. In reality, the diagonal diffusion coefficients may be negative,12 and sometimes the negative offdiagonal diffusion coefficients may be very large.8 The large negative diffusion coefficients may result in uphill diffusion, in which a species diffuses against its own concentration gradient, and multicomponent convection where the system is completely destabilized. These phenomena are not possible when the assumptions leading to (5.13) and (5.14) hold. To test the applicability of Cooper’s model and the common force model in a ternary system the mutual diffusion coefficients, self-diffusion coefficients, and thermodynamic activities must be all measured. Such measurements have only recently been made in a study of diffusion in ternary silicate melts.22 That study shows that in the system SiO2-Al2O3-CaO the mutual diffusion coefficients are better described by the common force model than by Cooper’s model.22 The specific negative crosscorrelation functions implicitly assumed in the common force model seem to be more appropriate for silicate liquids. Although the mutual diffusion coefficients of several other ternary solutions have been measured,12,38 tests of the empirical diffusion models for these systems are not possible because either self-diffusion coefficients or chemical activity data are unavailable at present. Although a recent simulation study attempted to calculate the mutual diffusion coefficients of a ternary silicate mixture,39 the results cannot be used for testing the diffusion models because the mutual diffusion coefficients were calculated incorrectly in that study. Finally, we point out that in our derivation the phenomenological coefficient matrix L is not constrained to be symmetric, though it is believed to be symmetric on the basis of Onsager’s theory.24 Equation 2.20 for ternary and higher-order systems provides a means of testing the Onsager reciprocal relations, LRβ ) LβR. From eqs 2.20 and 3.5 it can be seen that the two off-diagonal phenomenological coefficients in ternary mixtures, L12 and L21, will be equal only when
∫0∞ 〈Bj 1(t)‚Bj 2〉 dt ) ∫0∞ 〈Bj 2(t)‚Bj 1〉 dt
to test the Onsager reciprocal relation between thermal current and diffusive current in a binary simple system.42 Those simulation results, however, were inconclusive because of the difficulty in unambiguously defining the diffusive contributions to heat currents.43 In contrast, the evaluation of L12 and L21, as we proposed here, is conceptionally simple and straightforward. The resolution of this approach will be limited by the statistical precision with which the collective velocity correlation functions can be determined.
(5.15)
This clearly emphasizes the origin of the Onsager reciprocal relations in the principle of microscopic reversibility. Attempts to experimentally verify Onsager symmetry are subject to relatively large uncertainties.40,41 Efforts have also been made
Green-Kubo formulas for the mutual diffusion coefficients in n-component isothermal, isobaric systems have been derived. The diffusivity matrix is a product of a kinetic matrix and a thermodynamic matrix. The two matrices may be expressed in terms of time integrals of diffusion flux-correlation functions and spatial integrals of radial distribution functions, respectively. Results in the barycentric reference frame and volume-fixed reference frame are presented. Empirical diffusion models that link the mutual diffusion coefficients to the self-diffusion coefficients are discussed in relation to the statistical mechanics result. A zero net contribution from the cross-correlation functions yields the Darken or Hartley-Crank equation for binary systems (Cooper’s model for ternary systems), while a particular negative contribution from the cross-correlation functions results in Lasaga’s relation (or the common force model). The formulas for the phenomenological coefficients in ternary and higher-order systems provide a conceptually simple means to verify the validity of the Onsager reciprocal relations. Acknowledgment. We thank Yan Liang and Frank Richter for many stimulating discussions. This work was supported by NSF EAR-9304263 and by a gift from the Alcoa Foundation. Appendix A: Notation molar mass of component R partial molar volume of component R mole fraction of component R number of particles of component R chemical potential of component R per mole difference in specific chemical potential between components R and n mass density at equilibrium F m: F: number density at equilibrium N: total number of particles at equilibrium V: total volume T: equilibrium temperature P: equilibrium pressure self-diffusion coefficient of component R, eq 5.1 DR: radial distribution function of the pair of species R and β gRβ(r): spatial integral of gRβ(r), eq 3.10 ΓRβ: cR and cR: mass fraction or concentration of component R and its Fourier component at k ) 0 number density of component R and its Fourier component FR and at k ) 0, eqs 3.2a and 3.4a FR: BjR and BjR: particle current of component R and its Fourier component at k ) 0, eqs 3.2b and 3.4b JR: diffusion flux of component R and its Fourier component B JR and B at k ) 0, eq 3.5 D: mutual diffusivity matrix in the barycentric reference frame, eq 2.19 mutual diffusivity matrix in the volume-fixed reference DV: frame, eq 4.1 mR: V h R: xR: NR: µR: µRm:
Green-Kubo Formulas for Diffusion Coefficients T: L: C: C0: G: Ω:
J. Phys. Chem., Vol. 100, No. 13, 1996 5523
reference frame transformation matrix, eq 2.25 phenomenological coefficient matrix, eq 2.20 static concentration correlation function matrix, eq 2.14 C evaluated at k ) 0 thermodynamic matrix, related to the matrix C0 by eq 2.21 time-dependent concentration correlation function matrix, eq 2.13 matrix of quadratic particle number fluctuations, eq 3.8 matrix consisting of the integrals of radial distribution functions, eq 3.9 matrix of the derivatives of chemical potential with respect to particle number in the chemical ensemble, eq 3.12
N: B: A:
If we define the thermodynamic factor Q of a binary system as15
Q)
1 1 + x1x2F(Γ11 + Γ22 - 2Γ12)
(C1)
then its only independent chemical potential derivative is given by
( )
1 ∂µ1 kBT ∂x1
) T,P
Q x1
(C2)
Similarly, for ternary system the thermodynamic factors may be expressed in the following,
Appendix B: Expressions for 〈cr cβ〉 in Terms of the Integrals of Radial Distribution Functions Γrβ
Q11 )
Following the discussion in section 3.2, we can obtain microscopic expressions for the static concentration correlation functions. The general expression is:
1 [1 + x1x2F(Γ22 + Γ13 - Γ12 - Γ23) + x2x3F(Γ22 + Q′ Γ33 - 2Γ23)] (C3)
Q12 )
1 [x F(Γ + Γ13 - Γ12 - Γ23) - x3F(Γ33 + Γ12 Q′ 2 22 Γ13 - Γ23)] (C4)
Q21 )
1 [x F(Γ + Γ23 - Γ12 - Γ13) - x3F(Γ33 + Γ12 Q′ 1 11 Γ13 - Γ23)] (C5)
Q22 )
1 [1 + x1x2F(Γ11 + Γ23 - Γ12 - Γ13) + x2x3F(Γ11 + Q′ Γ33 - 2Γ13)] (C6)
R
〈c c 〉 ) β
mRmβxRN 2 Fm m j
xβm j (mR +
n
{xβ ∑ mγ2 xγ + m j 2δRβ γ)1
mβ) + xβF[m2nx2nΓnn n-1
-m j mnxn(Γβn + ΓRn) + n-1
j ∑ mγxγ(Γγβ + ΓγR) + m j 2ΓRβ + 2mnxn ∑ mγxγΓγn - m γ)1
γ)1 n-1 n-1
∑ ∑ mγmσxγxσΓγσ]}
(B1)
γ)1 σ)1
for R, β ) 1, 2, ..., n - 1. For n ) 2 the above expression is
〈c c 〉 ) 1 1
m21m22x1x2F2 N 4 Fm
[1 + x1x2F(Γ11 + Γ22 - 2Γ12)] (B2)
For n ) 3 the expression yields
〈c1c2〉 )
m21x1F2 N 4 Fm
{m22x2[x1 + x2 + x1x2F(Γ11 + Γ22 -
The denominator Q′ is given by
Q′ )
〈c c 〉 ) 2
m1m2x1x2F2 N 4 Fm
〈c2c2〉 )
m22x2F2 N 4 Fm
Following Kirkwood and Buff,27 we obtain the expressions for the derivatives of chemical potential with respect to mole fraction in terms of the radial distribution functions.
1 ∂µ1 kBT ∂x2
T,P,x1
∂µ2 kB T ∂x1
T,P,x2
∂µ2 kB T ∂x2 1
{m21x1 [x1 + x2 + x1x2F(Γ11 + Γ22 -
Appendix C: Expressions for (Dµr/Dxβ)T,P,xγ*β in Terms of the Integrals of Radial Distribution Functions Γrβ
) ) ) )
T,P,x2
1
[1 + x3F(Γ33 + Γ12 - Γ13 - Γ23)]} (B4)
( ( ( (
∂µ1 kB T ∂x1 1
{m1m2[x1 + x2 + x1x2F(Γ11 +
2Γ12)] + m23x3 [x2 + x3 + x2x3F(Γ22 + Γ33 - 2Γ23)] + 2m1m3x1x3 [1 + x2F(Γ22 + Γ13 - Γ12 - Γ23)]} (B5)
(C7)
where the matrix B is given by eq 3.9, and |B|Rβ is the cofactor of BRβ in the determinant |B|. Consequently, in terms of these thermodynamic factors, the chemical-potential derivatives with respect to mole fractions can be written as
Γ22 - 2Γ12)] + m1m3x3[1 + x1F(Γ11 + Γ23 - Γ12 - Γ13)] + m2m3x3[1 + x2F(Γ22 + Γ13 - Γ12 - Γ23)] m23x3
3
∑ ∑ xRxβ|B|Rβ x x x F R)1 β)1 1 2 3
2Γ12)] + m23x3[x1 + x3 + x1x3F(Γ11 + Γ33 - 2Γ13)] + 2m2m3x2x3[1 + x1F(Γ11 + Γ23 - Γ12 - Γ13)]} (B3) 1
3
1
Q11 x1
(C8)
) Q12
(C9)
) Q21
(C10)
Q22 x2
(C11)
)
) T,P,x1
Note that of the four thermodynamic factors Qij, only three are independent, and there exists the relation Q11 + (x2 + x3)Q21 ) Q22 + (x1 + x3)Q12. References and Notes (1) (2) 1203. (3) (4)
Kubo, R. J. Phys. Soc. Jpn. 1957, 12, 570. Kubo, R.; Yokota, M.; Nakajima, S. J. Phys. Soc. Jpn. 1957, 12, Zwanzig, R. J. Chem. Phys. 1964, 40, 2527. Zwanzig, R. Annu. ReV. Phys. Chem. 1965, 16, 67.
5524 J. Phys. Chem., Vol. 100, No. 13, 1996 (5) Rahman, A. Phys. ReV. B 1964, 6, A405. (6) Ciccotti, G.; Frenkel, D.; McDonald, I. R. Simulation of Liquids and Solids - Molecular Dynamics and Monte Carlo Methods in Statistical Mechanics; North-Holland: Amsterdam, 1987. (7) Oishi, Y. J. Chem. Phys. 1965, 43, 1611. (8) Liang, Y.; Richter, F. M.; Watson, B. E. Nature 1994, 369, 390. (9) Tyrrell, H. J. V.; Harris, K. R. Diffusion in Liquids; Butterworths: London, 1984. (10) Hofmann, A. W. In Physics of Magmatic Processes; Hargraves, R. B., Ed.; Princeton University Press: Princeton, 1980; pp 385-417. (11) Varshneya, A. K.; Cooper, A. R. J. Am. Ceram. Soc. 1972, 55, 418. (12) Vitagliano, V.; Sartorio, R.; Scala, S.; Spaduzzi, D. J. Solution Chem. 1978, 7, 605. (13) Watson, B. E.; Baker, D. R. In Physical Chemistry of Magmas; Kushiro, I., Perchuk, L., Eds.; Springer-Verlag: Berlin, 1991; pp 120151. (14) Cohen, C.; Sutherland, J. W. H.; Deutch, J. M. Phys. Chem. Liquids 1971, 2, 213. (15) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids; Academic Press: New York, 1986. (16) Darken, L. S. Trans. AIME 1948, 175, 184. (17) Bearman, R. J. J. Phys. Chem. 1961, 65, 1961. (18) McCall, D. W.; Douglass, D. C. J. Phys. Chem. 1967, 71, 987. (19) Cooper, A. R. Phys. Chem. Glasses 1965, 6, 55. (20) Lasaga, A. C. Geochim. Cosmochim. Acta 1979, 43, 455. (21) Richter, F. M. Geochim. Cosmochim. Acta 1993, 57, 2019.
Zhou and Miller (22) Liang, Y. Ph.D. thesis, University of Chicago, 1994. (23) DeGroot, S. R.; Mazur, P. Non-equilibrium Thermodynamics; Dover: New York, 1984. (24) Onsager, L. Ann. N. Y. Acad. Sci. 1945, 45, 241. (25) Kirkwood, J. G.; Baldwin, R. L.; Dunlop, P. J.; Gosting, L. J.; Kegeles, G. J. Chem. Phys. 1960, 33, 1505. (26) Mills, R.; Malhotra, R.; Woolf, L. A.; Miller, D. G. J. Phys. Chem. 1994, 98, 5565. (27) Kirkwood, J. G.; Buff, F. J. Chem. Phys. 1951, 19, 774. (28) Jolly, D. L.; Bearman, R. J. Mol. Phys. 1980, 41, 137. (29) Schoen, M.; Hoheisel, C. Mol. Phys. 1984, 33, 33. (30) Heyes, D. M. J. Chem. Phys. 1992, 96, 2217. (31) Trulla`s, J.; Padro´, J. A. Phys. ReV. E 1994, 50, 1162. (32) Douglass, D. C.; Frisch, H. L. J. Phys. Chem. 1969, 73, 3039. (33) Jacucci, G.; McDonald, I. R. Physics 1975, 80A, 607. (34) Zhou, Y.; Miller, G. H. Phys. ReV. E, in press. (35) Schoen, M.; Hoheisel, C. Mol. Phys. 1984, 52, 1029. (36) Hansen, J. P.; McDonald, I. R. Phys. ReV. A 1975, 11, 2111. (37) Gupta, P. K.; Cooper, A. R. Physica 1971, 54, 39. (38) Kim, H. J. Chem. Eng. Data 1982, 27, 255. (39) Kubicki, J. D.; Lasaga, A. C. Phys. Chem. Minerals 1993, 20, 255. (40) Miller, D. G. J. Phys. Chem. 1959, 63, 570. (41) Spera, F. J.; Trial, A. F. Science 1993, 259, 204. (42) MacGowan, D.; Evans, D. J. Phys. ReV. A 1986, 34, 2133. (43) Evans, D. J.; MacGowan, D. Phys. ReV. A 1987, 36, 948.
JP9533739
