10264
J. Phys. Chem. B 2008, 112, 10264–10271
Electrostatic Relaxation and Hydrodynamic Interactions for Self-Diffusion of Ions in Electrolyte Solutions J.-F. Dufrêche,*,† M. Jardat,† P. Turq,† and B. Bagchi*,‡ Laboratoire Liquides Ioniques et Interfaces Charge´es, case courrier 51, UniVersite´ P. et M. Curie - Paris 6, CNRS, 4 place Jussieu, 75252 Paris Cedex 05, France and Solid State and Structural Chemistry Unit, Indian Institute of Science, Bangalore, India 560012 ReceiVed: February 29, 2008; ReVised Manuscript ReceiVed: April 22, 2008
The concentration dependence of self-diffusion of ions in solutions at large concentrations has remained an interesting yet unsolved problem. Here we develop a self-consistent microscopic approach based on the ideas of mode-coupling theory. It allows us to calculate both contributions which influence the friction of a moving ion: the ion atmosphere relaxation and hydrodynamic interactions. The resulting theory provides an excellent agreement with known experimental results over a wide concentration range. Interestingly, the mode-coupling self-consistent calculation of friction reveal a nonlinear coupling between the hydrodynamic interactions and the ion atmosphere relaxation which enhances ion diffusion by reducing friction, particularly at intermediate ion concentrations. This rather striking result has its origin in the similar time scales of the relaxation of the ion atmosphere relaxation and the hydrodynamic term, which are essentially given by the Debye relaxation time. The results are also in agreement with computer simulations, with and without hydrodynamic interactions. 1. Introduction Electrolyte solutions exhibit unusual concentration dependences of ionic conductivity, self-diffusion of ions, and viscosity of the solution.1 The well-known Debye-Hu¨ckel-Onsager law of ionic conductivity2 and the Onsager-Fuoss-Falkenhagen law of viscosity3 both provide a square-root concentration dependence in the limit of low concentrations.4 Although derivation of these laws constitute one of the significant intellectual triumphs of 20th century, unfortunately they are valid only at very low concentration, often at less than 0.05 M. However, electrolyte solution shows a remarkable variety of concentration dependences which has till now eluded a fully microscopic explanation and a quantitative understanding. At very low concentration where the distance between ions is large, transport properties are determined by the ion atmosphere relaxation and electrophoretic effects.5 Here the Debye length6 is much larger than the size of the solvent molecules and the microscopic aspects of solvent effects can be safely neglected. Difficulty appears when the concentration becomes large such that the Debye length is no longer much larger than the solvent molecular size. Strictly speaking, in this limit, one should develop a theory which takes into account all the degrees of freedom of solvent molecules at the same level as ions.7 However, since we are mostly interested in aqueous solutions, there is a formidable difficulty in developing such a theory, the solvent itself having complex structural and dynamical properties. In the absence of such a fully microscopic theory, attempts have been made8–15 to develop an effective theory where static intermolecular correlations between ions are taken into account with the help of a statistical mechanical model, like the mean spherical approximation (MSA).16,17 The dynamical correlations of the solvent (still treated as a continuum) have been included * To whom correspondence should be addressed. E-mail:
[email protected] and
[email protected]. † Universite ´ P. et M. Curie. ‡ Indian Institute of Science.
in a most recent theory by using a mode-coupling theory which allows a self-consistent treatment of diffusion and conductivity.18–21 The above approach has met with considerable success in describing the concentration dependence of ionic conductivity and self-diffusion of ions at large concentration, even up to 2 M solutions in some cases. The above theoretical treatment, however, ignores a potentially important effect, which is hydrodynamic interaction between ions. The hydrodynamic interactions (HIs) are a purely kinetic effect propagated through the solvent molecules and known to be of importance in the dynamics of colloids.22 In the case of ions in electrolyte solutions, HIs can be of significance because pair interaction is also of long range and the ion atmosphere relaxation can be quite slow at intermediate concentration. Thus, there should be a window of concentration neither too low nor too high where HIs can play an important role. Recent computer simulation studies of self-diffusion in electrolyte solution have indeed confirmed such expectation, as discussed later.23 The objective of this work is to develop a mode-coupling theory of the hydrodynamic interaction term which is projected into the density relaxation of the electrolyte, the latter being given by the dynamic structure factor. The reason why such a theory can be meaningfully developed is that the dyadic tensor that enters in the Oseen tensor description of HIs essentially contains two density terms, as eloquently shown by Denny and Reichmann in their treatment of third- and fifth-order nonresonant Raman spectra in liquids.24,25 As it will be shown in the following, the resulting modecoupling theory of self-diffusion contains three terms: the first one for the ion atmosphere relaxation, the second one due to hydrodynamic interaction, while a third term comes from the cross-correlation between the first two. It turns out that the third term plays a highly interesting role in reducing electrolyte friction and enhancing diffusion, as demonstrated in more detail later. The resulting theory is found to be in agreement both with known experimental results and with computer simulations
10.1021/jp801796g CCC: $40.75 2008 American Chemical Society Published on Web 07/08/2008
Self-Diffusion of Ions in Electrolyte Solutions
J. Phys. Chem. B, Vol. 112, No. 33, 2008 10265
(Brownian dynamics, BD) which have been performed specifically to test the amplitude of HIs in ionic self-diffusion. The organization of the rest of the paper is as follows. In section 2, we describe the model and define many useful quantities which are used in section 3, where we present the theoretical approach. The ion atmosphere relaxation and HIs terms are calculated in section 4. The Brownian dynamics simulations are described in section 5. The numerical results are presented in section 6, and our conclusions are summarized in section 7. 2. Basic Model and Definitions Let us first consider an electrolyte solution which is composed of three species: the solvent, cation, and anion. The dynamics of the solute particles will be studied at the limit where the solvent can be considered as a continuum of dielectric constant εr and viscosity η. Within this continuous solvent description the ions interact through a solvent-averaged effective potential (McMillan-Mayer description26,27) which is a pair potential for dilute solutions. It is the sum of two terms SR VRβ(r) ) VRβ (r) +
ZRZβe2 4πε0εrr
(1)
where r is the distance between two particles of species R and β. For the sake of clarity, throughout this paper the Greek letters (R, β, γ, δ) will refer to the species of particles, whereas the Latin letters (i, j) will denote the particles themselves. e is the elementary charge, and ZR is the charge of the ion species R. SR The short-range interaction potential VRβ (r) (which can be 28,29 obtained from all-atom simulations ) will be modeled by a hard-sphere repulsion. This so-called primitive model3 is very commonly used in the studies of structure30 and dynamics14 of electrolyte solutions. The Fourier transform of the time-dependent number density of species R FR(r, t), with r being the position at time t, reads
FR(k, t) )
∫-∞+∞ FR(r, t)eikr dr
(2)
The intermediate scattering function FRβ(k, t) is given by31
FRβ(k, t) )
1
√NRNβ
〈δFk(t)δF-k(0)〉
(3)
where δFk(t) is the Fourier transform of FR(r, t) - FR, the density fluctuation of species R. The Laplace transform of FRβ(k, t) is
FRβ(k, z) )
∫0+∞ FRβ(k, t)e-zt dt
(4)
3. Hydrodynamic Interactions and Relaxation Effect 3.1. Significance of Hydrodynamic Interactions. Since the pioneering works of Debye, Hu¨ckel, and Onsager (DHO), it is well known that the hydrodynamic interactions (HIs) together with (electrostatic) relaxation are one of the two effects which control the ion transport properties in electrolyte solutions.2,3,5,32 For example, in the case of the conductivity,12,33 the HIs yield roughly one-half of the variation of this transport coefficient as a function of the concentration for simple electrolyte solutions. In the case of ion self-diffusion though34 they do not influence the variation of that transport coefficient as a function of the concentration for the regime of high dilution, which corresponds to the DHO limiting laws, where the variation is proportional to the square root of the concentration. It arises from the fact that self-diffusion corresponds to an experiment where there is only a drag force on the tagged tracer particles. In fact, there is
no net force on the surrounding particles of the tagged particles (because the latter are infinitely diluted) and the tracer particles are not submitted to hydrodynamic interactions. However, this description of the phenomenon is only valid in the linear regime, where the transport coefficient is proportional to the square root of the concentration. Even if for higher concentration HIs themselves are found to be negligible,21 they are evidence from Brownian dynamics simulations23,35 that they can still influence quantitatively the ion self-diffusion coefficient at high concentration. This nonlinear effect arises from the fact that the (electrostatic) relaxation is modified by HIs. Indeed, when the tracer particle moves, the solvent molecules around it influence the velocities of the surrounding ions of the ionic atmosphere, so that the total relaxation force is modified. The theory developed here makes the origin of the implicit role of HIs clear. 3.2. Electrolyte Friction. The microscopic modeling of such a nonlinear effect is far from being obvious. There are several nontrivial aspects of this problem. First, we need to translate the ion-ion interaction term into a time correlation formalism. Second, we need to describe the interplay between the HIs and the electrostatic relaxation term. Third, we need to develop a self-consistent theory. Mode-coupling theory (MCT) shows that the microscopic friction ζs(z) can be expressed as18,36–38
ζs(z) ) ζSR(z) + δζmic(z)
(5)
The first term ζSR(z) corresponds to short-range interactions, and δζmic(z) corresponds to long-range interactions which depend on the electrolyte concentration. ζSR(z) is the sum of several effects
ζSR(z) ) ζb(z) + ζFF(z) + ζF,DF(z)
(6)
where ζb(z) denotes the binary collision term, ζFF(z) represents the isotropic density fluctuations39 (which are present even in a nonpolar solvent), and ζF, DF(z) denotes the dielectric friction of the solvent.40,41 These quantities could be, at least in principle, calculated in a fully microscopic theory. They express the ion-solvent interactions. Consequently, they are independent of the electrolyte concentration for moderately concentrated solutions (as long as the hydration is not significantly modified). Thus, they simply represent the friction at infinite dilution when only ion/solvent interactions are present. At the Brownian limit for the ions (which is a very satisfying approximation for electrolyte solutions even for small ions because of the hydration5,13,42), their relaxation time is much smaller than the relaxation time of the ion degrees of freedom. The corresponding velocity autocorrelation function whose Laplace transform is
Z0s (z) )
kBT msz + ζSR(z)
(7)
has a time integral which corresponds to the self-diffusion coefficient of the ion at infinite dilution Ds0
D0s ) Z0s (z ) 0) )
kBT ζSR(z ) 0)
(8)
In the second term of eq 5, δζmic(z), is the electrolyte friction which depends on the electrolyte concentration c. It is equal to zero at infinite dilution and much smaller than msz + ζSR(z) for dilute solutions. Its relaxation time is the one of the ionic atmosphere and can be roughly identified to the electrostatic Debye relaxation time τD,21,43,44 which is typically on the order of 50 ps for a simple 1-1 aqueous electrolyte at 1 M. In fact,
10266 J. Phys. Chem. B, Vol. 112, No. 33, 2008
Dufrêche et al.
its relaxation time is much higher than the relaxation time of ζSR(z). We have
Z(z) )
kBT msz + ζ(z)
kBT kBT δζ (z) ≈ msz + ζSR(z) (m z + ζ (z))2 mic s SR ≈
where ζ0 ) ζ0(z ) 0). 3.3. Linear Response Theory at the Smoluchowski Level of Description. The microscopic electrolyte friction δζmic(z) can be calculated using the Kirkwood formula31,36
1 3kBT
ions -zt dt ∫0+∞ 〈Fions s (t)Fs (0) 〉 e
(10)
where Fions s (t) is the time-dependent force exerted on the tagged ion due to its interactions with all other ions in the solution. Then the total velocity autocorrelation function reads
( )
0 2 1 Ds ions Z(t) ) Z (t) 〈Fions s (t)Fs (0)〉 3 kBT 0
(11)
The first term stands for the solvent/ion interactions and decays much faster than the second one. When the hydrodynamic interactions are neglected, we already proposed a method to calculate the second one, which represents the relaxation effect.21 Here the issue consists of evaluating this term by taking into account the HIs. This expression (eq 11) of Z(t) is consistent with the Smoluchowski theory of electrolyte solutions,11,14,33,42,45 where the Kubo-like formula46,47 of the self-diffusion coefficient when no HIs are taken into account reads23,48
Ds ) D0s -
1 3
∫0
+∞
( )
D0s 2 MM 〈Fs (t)FMM s (0)〉 kBT
dt
(12)
with FsMM(t) being the McMillan-Mayer force on the tracer particle averaged over the configurations of the solvent. Then eqs 11 and 12 are the same when the HIs are neglected if Fsions can be identified to FsMM. At first sight, when hydrodynamic interactions are taken into account, the two formulas do not coincide anymore. Indeed, the generalized expression of eq 12 with HIs is23
Ds ) D0s -
1 3
∫0+∞ dt〈Us(t)Us(0)〉
(13)
Us(t) is the hydrodynamic velocity of the tracer, which is the velocity of the tracer that can be calculated from the hydrodynamic equations, i.e.
Us(t) ) βD0s FMM + USs s
(14)
where β ) 1/kBT and
∑ TsiFiMM
{
(9)
kBT kBT - 0 2 δζmic(z) msz + ζSR(z) (ζ )
δζmic(z) )
are supposed to be point-like from the hydrodynamic point of view, so, consistent with the previous study of transport of j si is the Oseen tensor electrolyte solutions, we will consider that T defined by
(
)
1 rXr II + 2 8πηr r 1 (r · f)r j f+ 2 i.e. Tsif ) 8πηr r
j ) T si
(
)
(16)
where r ) rs - ri is the distance vector between the particles and η the viscosity of the solvent. Equation 13 is more explicit than eq 11 because it explicitly involves the hydrodynamic interactions. Thus, this Smoluchowski-level formula will be our basic formula for calculation of the self-diffusion coefficients. 3.4. Velocity Correlation Function. In this section we are going to show from a deeper analysis that eq 13 is actually equivalent to eq 11 if we correctly write down the force Fsions(t) exerted on the tagged ion due to its interactions with all other MM ions. Fions because of the hydrodynamic s (t) is different from Fs interactions. There is an alternative method which can be used to obtain eq 13. Instead of using the Smoluchowski equation, it is possible to start from the Langevin equations, which describe the motion of the solute particles. When no HIs are considered, we have
m
dvs ) -ζ0vs + FiMM(t) + fs(t) dt
(17)
where vs is the velocity of the tracer particle and fs(t) is the random force which represents the collision with solvent molecules. It evolves much faster than the other forces, so that it can be modeled by a Gaussian white noise which satisfies the fluctuation-dissipation theorem 〈fs(t) · fs(0)〉 ) 6kBTζ0δ(t). When HIs are taken into account, the friction of the solvent -ζ0vs is modified by the presence of solute particles. A simple treatment consists of assuming that vs has to be replaced by vs minus the velocity of the solvent modified by the surrounding particles UsS. Thus, the Langevin equation becomes
m
dvs ) -ζ0(vs - USs (t)) + FMM s (t) + fs(t) dt
(18)
with the Langevin eqution, eq 18, being linear, it is possible to calculate vs(t) as a function of FsMM(t) and fs(t). Then, by averaging over the initial position of the tracer particle we get the velocity correlation function. More precisely straightforward calculations shows that
Z(t) )
kBT ζ0t⁄m 1 〈(ζ0USs (t) + FMM e s (t)) 02 m 3ζ × (ζ0USs (0) + FMM s (0))〉 (19)
Since 1/ζ0 ) Ds0/kBT and Ds ) ∫+∞ 0 Z(t) dt, we finally get
Ds ) D0s -
D0s 3kBT
∫0∞ dt〈(ζ0USs (t) + FMM s (t))
(15)
× (ζ0USs (0) + FMM s (0))〉 (20)
represents the hydrodynamic velocity of the solvent which arises from the motion of the surrounding solute particles at the place where there is the tracer particle. i denotes the different solute particles (but the tracer). FMM is the McMillan-Mayer force i j si is the hydrodynamic tensor between s and on particle i and T i which depends on the position of the solute particles. The ions
which is exactly equivalent to eq 13. Thus, hydrodynamic interactions modify the relaxation term because we have to consider also the relaxation of the solvent flux which is created by the relaxation of the ionic atmosphere. We will use eq 13 as the starting point of our treatment. It can be derived from eq 11, obtained from the microscopic friction, if the hydrodynamic force ζ0UsS(0) is taken into account in the time-dependent force
USs )
i
Self-Diffusion of Ions in Electrolyte Solutions
J. Phys. Chem. B, Vol. 112, No. 33, 2008 10267
Fsions(t) exerted on the tagged ion due to its interactions with all other ions MM 0 S Fions s (t) ) Fi (t) + ζ Us (t)
(21)
The latter expression can be seen as a direct consequence of the Langevin equation with HIs (eq 18) which separates the effects of the ions and the effects of the solvent. 3.5. Coupling between Hydrodynamics Interactions and Electrostatic Relaxation Effect. The general expression of the self-diffusion coefficient of an ion eq 13 consists of several terms which are as follows
Ds )D0s -
( )∫
0 1 Ds 3 kBT
2
+∞
0
simple product of 2-body correlation functions, we obtain an expression for Crr(t). The latter reads in the Fourier space
Crr(t) )
kB2T2
∑ √FRFβ ∫ dkFs(k, t)csR(k)csβ(k)FR,β(k, t)k2
(2π)3 R,β
(25) The calculation of this term has been discussed before.21 4.2. HIs Relaxation Coupling Term. The third term of eq 22 can be calculated by a similar method
dt〈Fs(t)Fs(0)〉
C1hr(t) )
∫+∞ dt〈∑ TsiFi(t)Fs(0) 〉 3kBT 0
-
D0s 3kBT 1 3
∫0
+∞
〈
dt Fs(t)
∑
TsiFi(0)
i
)
1 V
∫ dr1 ∫ dr2 ∫ dr3∑ 〈Fs(r2, t) ×
(22)
In this last formula, all the forces denoted by F correspond to McMillan-Mayer forces, although we removed for simplicity all the MM superscript on the forces. The first term is the selfdiffusion coefficient at infinite dilution which represents the fast decaying mode associated with the solvent. The second term is the simple relaxation term which represents relaxation of the electrostatic force. The second line stands for the coupling between the relaxation of the electrostatic force and the hydrodynamic interactions. The last term represents the relaxation of the hydrodynamic interactions. Consequently, the present treatment of the influence of the HIs on the ion selfdiffusion coefficient shows that the HIs not only modify the electrostatic relaxation around the tracer particle but also influence Ds directly. Thus, eq 22 contains four different kinds of contributions. Interestingly, relaxation of the hydrodynamic velocity of the solvent couples with the electrostatic relaxation to give rise to a cross-term which we show to be significant. 4. Calculation of the Ion Atmosphere Contribution Modified by HIs 4.1. Pure Relaxation Term. The second term of eq 22 can be calculated from time-dependent density functional theory.18,21 Indeed
Crr(t) ) 〈Fs(t) · Fs(0) 〉 )
∫ dr1 ∫ dr2〈Fs(r1, t)Fs(r2, 0)〉 (23)
where Fs(r, t) is the force exerted to the ion per unit of volume. An expression for Fs(r, t) can be obtained from time-dependent density functional theory7
Fs(r, t) ) kBTFs(r, t)gradr
(26)
γ
Substitution of eq 24 into eq 26 yields
j
i
s
jT (r - r )F (r , t)F (r , 0)〉 sγ 2 3 γ 3 s 1
〉
∫0+∞ dt〈∑ TsiFi(t)∑ TsjFj(0)〉
si i
i
D0s
i
-
〈∑ T F (t)F (0)〉
∑ ∫ dr ′ csR(r, r ′ )δFR(r ′ , t) R
(24) where δFR(r, t) ) FR(r, t) - FR is the density fluctuation of particles of species R and csR(r, r′) is the direct correlation function between species R and the tagged ion s. By introducing eq 24 into eq 23 and using the Gaussian approximation which consists of approximating the N-body correlation function by a
C1hr(t) )
(kBT)2 V
∫ dr1 ∫ dr2 ∫ dr3 ∫ dr4 ∫ dr5 ∑
R,β,γ
〈Fs(r2, t)Tsγ(r2 - r3)Fγ(r3, t)gradr3cγβ(r3, r4)δFβ(r4, t) × Fs(r1, 0)gradr1csR(r1, r5)δFR(r5, 0)〉 (kBT)2 ≈ V
∫ dr1 ∫ dr2 ∫ dr3 ∫ dr4 ∫ dr5 ∑
R,β,γ
j (r - r )F (r , t)F (r , 0)〉 〈Fs(r2, t)T sγ 2 3 γ 3 s 1 × 〈 gradr3cγβ(r3, r4)δFβ(r4, t)gradr1csR(r1, r5)δFR(r5, 0)〉 (27) where the inner part of the angular bracket has been approximately factorized. Then we obtain
C1hr(t) ≈
(kBT)2 V
∫ dr1 ∫ dr2 ∫ dr3 ∫ dr4 ∫ dr5 ∑
R,β,γ
Fγωsγ(r2 - r3)Gss(r2 - r1, t) × 〈gradr3cγβ(r3, r4)δFβ(r4, t)gradr1csR(r1, r5)δFR(r5, 0)〉 (28) where Gss(r, t) is the self Van Hove function and ω j sγ(r) ) j sγ(r) with gsγ(r) the pair correlation function between gsγ(r))T s and γ. We now substitute eq 24 into eq 28, and the Fourier transform yields
C1hr(t) ) kB2T2
∑ √FRFβFγ ∫ dkFs(k, t)csR(k)cγβ(k)FR,β(k, t)kωsγ(k)k
(2π)3 R,β,γ
(29) 2 Because of the time reversal symmetry the next term Chr (t) ) 1 j 〈Fs(t)ΣiTsiFi(0)〉 is equal to Chr(t), so that we can combine the two
10268 J. Phys. Chem. B, Vol. 112, No. 33, 2008
Dufrêche et al.
〈∑ T F (t)F (0) 〉 + 〈F (t)∑ T F (0)〉
Chr(t) )
si i
s
s
i
i
2kB2T2 3
)
FRβ(k, z) ) [-z + DR(z)k2]-1SRβ(k)
si i
(2π)
+
∑ √FRFβFγ ∫ dkFs(k, t)csR(k)
R,β,γ
× cγβ(k)FRβ(k, t)k · ωsγ(k)k (30) 4.3. HIs Relaxation Term. The last term of eq 22 reads
Chh(t) )
〈∑ T F (t) · T F (0)〉 si i
(31)
sj j
i,j
After a similar calculation one obtains
kB2T2
Chh(t) )
∑ √FRFβFγFδ ∫ dkFs(k, t)
× cδR(k)cγβ(k)FRβ(k, t)ωsγ(k)k · ωsδ(k)k (32) 4.4. Numerical Calculations: Practical Aspects. The selfdiffusion coefficient and velocity autocorrelation function can be obtained from eqs 25, 30, and 32 since we have
( )
0 2 D0s 1 Ds 1 Z(t) ) Z (t) C (t) C (t) - Chh(t) (33) 3 kBT rr 3kBT hr 3
We next discuss the terms that are required for the practical calculation Ds. First we have to calculate the vertex term
ωsγ(k)k )
∫ dr gsγ(r)Tsγ(r)e
ikr
k
(34)
We model the HIs by the Oseen tensor. We can remove the part which corresponds to far-field hydrodynamic field.48–51 Indeed
ωsγ(k)k )
∫ dr hsγ(r)Tsγ(r)eikrk + ∫ drTsγ(r)eikrk
(35)
and the second term which is the Fourier transform of the Oseen tensor (eq 16) applied to k gives a zero contribution
(
)
(
)
1 rXr 1 kXk II + 2 eikrk ) 2 II - 2 k ) 0 ∫ dr 8πηr r ηk k
(36)
We finally obtain
ωsγ(k)k )
k η
∫0+∞ rhsγ(r) kr cos(kr) + k(k3rr3
2 2
- 1)sin(kr)
dr (37)
Second, we have to know the self-dynamic structure factor of the tagged ion. Consistent with mode-coupling theory, we considered the hydrodynamic limit of this quantity
Fs(k, t) ) e-Dsk t 2
2
∑ √FRFγCRγ(k)Fγβ(k, z) (39)
-z + DR(z)k2 γ)1
where the frequency-dependent diffusion coefficient of species R DR(z) is approximated by the self-diffusion coefficient calculated by the present theory. This expression (eq 39) neglects the effect of the distinct diffusion coefficients, 10,11 which are small for dilute solutions. SRβ(k) ) FRβ(k, t ) 0), where SRβ(k) is the partial static structure factor between species R and β. SRβ(k) is related to the Fourier transform of the pair correlation function hRβ(k) by the following relation
SRβ(k) ) δRβ + √FRFβhRβ(k)
(2π)3 R,β,γ,δ
0
DR(z)k2
(38)
It is worth noting that this quantity which is required to calculate Ds depends itself on Ds. Thus, we have to make a selfconsistency loop to obtain the correct result: calculation of the self-diffusion coefficients Ds has to be performed several times by reinjecting the previous calculated values at every stage until the result converges. Third, we need to know the ionic van Hove function FRβ(k, t). Use of time-dependent density functional theory leads to the following equation for the frequencydependent van Hove function 18,52
(40)
The four coupled equations (R, β ) 1, 2) as given by eq 39 can be solved analytically to obtain the frequency dependence of the ionic van Hove functions. Finally, we require the solutions of the static structure factors and the direct correlation functions. The direct correlation functions are related to the static structure factors by the Ornstein-Zernike equations.31 The solutions of the pair correlation functions required for the calculation of these quantities have been obtained from the solution of the mean spherical approximation (MSA)16,17
hRβ(k) ) -
4π ˜ I G (ik) k m Rβ
(41)
˜ Rβ is the Laplace transform used in the Baxter-Wertheim where G factorization method for charged systems.17 We considered the whole analytical solution of the MSA with different sizes so that the Stillinger-Lovett sum rule is exactly verified. The integrals over wavevectors have been calculated numerically up to several times the inverse of the smaller size of the ions. 5. Brownian Dynamics Simulations On the mesoscopic time scale of Brownian dynamics, the motions are described in the position space only by the following stochastic equation for the displacement ∆r of N particles included in the simulation box from time t to t + ∆t53
(
r(t + ∆t) - r(t) ) βD · F +
∂ · D ∆t + R ∂r
)
(42)
where β ) 1/kBT. Here, the particles are supposed to be spherical, without rotational degree of freedom; ∆t is the time increment, r is the 3N-dimensional configuration vector, and F describes the forces acting on the particles at the beginning of the step. R is a random displacement, chosen from a Gaussian distribution with zero mean, 〈R〉 ) 0, and variance 〈RRT〉 ) 2D∆t. Hydrodynamic interactions (HIs) between particles are introduced via the configuration-dependent 3N × 3N diffusion matrix D. When HIs are neglected, D becomes diagonal and constant with eigenvalues equal to the self-diffusion coefficients of ions at infinite dilution. The Rotne-Prager tensor is used to calculate the diffusion matrix with some approximations;23,54 notice that the divergence of the diffusion matrix is strictly zero in that case. The Rotne-Prager tensor is used instead of the Oseen tensor because the diffusion matrix has to be positive-definite to compute the random displacements with an adequate distribution. The matrix obtained using the Oseen tensor is not always definite-positive conversely to the Rotne-Prager one. The first approximation consists of assuming that HIs are pairwise additive. This is reasonable for the systems we studied because they have a low packing fraction (less than 3%). Second, we consider hydro-
Self-Diffusion of Ions in Electrolyte Solutions dynamic interactions between particles included in the simulation box only. This constitutes a strong approximation. This approximation has been a posteriori justified several times for ionic solutions as we obtained very similar transport coefficients from two simulations which differ only from the total number of particles (see, for example, ref 55). The parameters of the simulations are the following. Aqueous solutions of KCl and NaCl are simulated at 298.15 K for three different concentrations: 0.25, 0.64, and 1.0 mol L-1. In each simulation, 300 ions are placed in a cubic box with periodic boundary conditions. To compute the soft core interactions we use a spherical cutoff of one-half a box length, applying the minimum image convention. Coulomb interactions are computed using the Ewald summation technique with the conducting boundary condition. The simulation time step is in the range 0.07-0.1 ps, depending on the electrolyte concentration. Dynamical properties are calculated by averaging over five trajectories, each containing 8.2 × 104 time steps, while the equilibration runs contain at least 2 × 105 time frames. In order to calculate the averaged dynamical properties over several independent trajectories, we perform one long run (about 2 × 105 time steps) without HIs and without saving any position or velocity between two production runs. Coulomb interactions are computed using the Ewald summation technique with the conducting boundary condition. There are a few differences between the model used for BD calculations and that used for MCT. First, the level of description for HIs is not exactly the same (Rotne-Prager vs Oseen). Second, the BD algorithm requires continuous forces so that continuous McMillan-Mayer potentials (i.e., no hard-sphere repulsion) are necessary. The definition of the MM potentials for BD follows the procedure described in ref 48: soft potentials of BD have been chosen so that they yield the same osmotic coefficients as hard potentials of MCT based on the MSA.
J. Phys. Chem. B, Vol. 112, No. 33, 2008 10269
Figure 1. Self-diffusion coefficients of aqueous electrolyte solutions at 25 °C as functions of the square root of concentration: (solid lines) theory with HIs; (dashed lines) theory without HIs; (/) experimental values; (∆) BD with HIs; (∇) BD without HIs.
6. Results and Discussion The theory has been applied to simple aqueous electrolyte solutions at 25 °C: εr ) 78.3 and η ) 0.89 × 10-3 Pa s. The only unknown parameters are the (hydrated) diameters of the ions. They have been taken from refs 13 and 14. These diameters have been obtained from the investigation of the equilibrium properties, but they are also able to reproduce the various transport properties self-consistently. Therefore, there is no adjustable parameters in our treatment and the results are real predictions. In Figure 1, the self-diffusion coefficients of the ions are plotted as functions of the square root of the concentration for NaCl, KCl, and LiCl solutions. The diameter of the ions are σNa+ ) 3.05 Å, σK+ ) 2.95 Å, σLi+ ) 4.35 Å, and σCl) 3.62 Å. The experimental data have been taken from the critical analysis given in ref 56. For NaCl, the theoretical predictions are in excellent agreement with the experimental results even at high concentrations. For LiCl and KCl, the agreement is also good, even if for KCl the theory slightly underestimates the self-diffusion coefficients for the highest concentrations. It is worth noting that it is possible to recover the experimental data by increasing the size of the ions. The comparison between the curves with and without HIs shows that HIs increase the self-diffusion coefficients. The agreement between the theoretical predictions and the experimental values is better when HIs are taken into account, as it can be shown for KCl. The comparison between MCT and BD is also interesting. When HIs are neglected, the agreement is excellent. When they are taken into account, BD confirms that HIs enhance diffusion.
Figure 2. Velocity autocorrelation function of Cl - in NaCl aqueous electrolyte solutions at 25 °C as a function of the time (SI unit) for two electrolyte concentrations: (solid lines) theory with HIs; (dashed lines) theory without HIs terms; (∆) BD with HIs; (∇) BD without HIs.
The effect appears to be a little more important for BD, but the difference is close to the error bar of the simulation results. The difference seems to be systematic though, and it may be due to the differences between the interaction potentials (soft core for BD and hard core for MCT) and between the treatments of the HIs (Rotne-Prager tensor for BD and Oseen tensor for MCT). We recover the same agreement between MCT and BD for the velocity autocorrelation function (Figure 2). HIs reduce the magnitude of the relaxation force but do not modify the relaxation time. The self-consistency is important for MCT. For example, in our case it modifies the magnitude of the effects by as much as a factor of 10% at high concentration. Furthermore, it expresses the coupling of the different effects, e.g., the fact that HIs modify the electrostatic relaxation. The influence of the different terms as a function of concentration is analyzed in Figure 3. The electrostatic relaxation
10270 J. Phys. Chem. B, Vol. 112, No. 33, 2008
Dufrêche et al.
Figure 3. Self-diffusion coefficients Na + in NaCl aqueous electrolyte solutions at 25 °C as a function of the square root of concentration. Influence of the various terms on the final result.
effect decreases self-diffusion, but its influence is less important at high concentration (the curve is above the limiting law) because of the size of the ions whose repulsion compensates for a little of the electrostatic attraction. The figure shows that the predominant term for HIs is the coupling term, which is positive. Indeed, the motion of the tracer particle moves hydrodynamically the ions of its ionic atmosphere so that they go to the same direction, which reduces the electrostatic force. The effect is significant because the long-range of HIs is coupled to the long range of the electrostatic relaxation forces. On the other hand, the pure hydrodynamic relaxation term decreases the self-diffusion, but this effect is negligible. It is worth noting that the pure electrostatic curve is slightly different than the curve where there is no HIs. This is an effect of the selfconsistency. The electrostatic relaxation effect is modified by the HIs: they reduce its amplitude because they increase the self-diffusion coefficient. Nevertheless, this effect is almost negligible. The main reason why HIs enhance self-diffusion is that there is a coupling between hydrodynamics and electrostatic forces. This enhancement is better understood if one combines the various terms of the total velocity correlation (eq 33) as follows
Z(t) ) Z0(t) kB2T2
∑ √FRFβ ∫ dkFs(k, t)FRβ(k, t)µsR(k)k · µsβ(k)k (43)
3(2π)3 Rβ
where the generalized mobility µ j
µsR(k)k )
(
D0s c (k) + kBT sR
sR(k)k
is
)
∑ FγcγR(k)ωsγ(k) k γ
(44)
Thus, in the long-time part of the velocity correlation function the vertex terms ω j (k) cancel out a part of the relaxation effect. Microscopically it corresponds to the following mechanism. When an external force is applied on a tracer ion, its ionic atmosphere is not spherical anymore and its deformation produces the electrostatic relaxation effect. However, at the same time, the motion of the tracer stirs the solvent molecules which drag the surrounding ions. In fact, when HIs are taken into account the deformation of the ionic atmosphere is reduced and the relaxation effect is lower. This analysis is confirmed if we look at the velocity correlation functions themselves (Figure 4). The characteristic time of the three correlation functions Crr(t), Chr(t), and Chh(t)
Figure 4. Velocity autocorrelation function of Na + in NaCl aqueous electrolyte solutions at 25 °C as a function of the time (SI unit): (solid line) full terms but the (fast) short-ranged term Z0(t) (which corresponds to Ds0), (dotted line) electrostatic relaxation term; (dashed line) pure hydrodynamic relaxation term; (dotted-dashed line) electrostatic-hydrodynamic coupling term.
is roughly the same; it corresponds to the Debye time. Thus, different plots at different concentration (not shown) indicate that the higher the concentration, the faster is the long time decay. As we already indicated in another paper,21 these results can explain the observed discrepancy between the self-diffusion coefficients measured by time-of-flight neutron scattering and NMR or tracer methods.57 The typical time scale of neutron scattering experiments is typically 20 ps, which is less than the relaxation time of the ionic atmosphere. Consequently, the values obtained do not account for the whole relaxation effect. They are found to be greater than those measured by long time methods such as NMR. However, they do depend on the concentration because part of the electrostatic relaxation is taken into account. The agreement with the experimental values is good, yet the model of the short-range part of the potential of mean force we used is very simple since the ions are only characterized by their hydrated diameters. It is very different from the effective potential obtained from explicit solvent calculations.28,29 In particular, in molecular descriptions, large distributions of opposite charge ions are often found at the contact distance and the force exhibits large oscillations which are reminiscent of the layering of the solvent between the ions. Our work suggests that such molecular effects do not appear to play an very important role for self-diffusion in the case of dilute solutions. Several explanations can be proposed. First, the fact that large cation-anion distributions can be found at the contact distance indicates that pairs of ions exist in the solution. This association effect can be taken into account in electrolyte theory.58 However, as long as the concentration of the pairs is low compared to the total ion concentration, they can be neglected. Second, selfdiffusion in liquids is strongly coupled to the hydrodynamical modes (that is the essential idea of mode-coupling theory), but short time collisions represent only a small part of the correlation function.36 Third, for dilute electrolytes, self-diffusion is mainly driven by the long-range part of the potential (Onsager limiting law), so that the short-range part can be seen as a perturbation effect. The situation could be totally different close to an interface where the solvent hydration and the dehydration of the ions can play an important role. However, in the case of bulk electrolyte solutions, it is not to bad an approximation to model the short-range part as a hard-sphere repulsion, as it is commonly done in perturbation theory of fluids. 7. Conclusion Let us first summarize the main results of this paper. We developed a mode-coupling theory of self-diffusion that includes
Self-Diffusion of Ions in Electrolyte Solutions contributions to the friction from both the interion pair interactions and the hydrodynamic interactions due to the solvent dynamic effects. Interactions among the ions are included through a time-dependent density field. The resulting modecoupling theory was described earlier. The effects of hydrodynamic interactions are taken into account through the wellknown Oseen tensor. A major component of this work is the mode-coupling theory expression of the force term (given by the Oseen tensor) that we derived by recognizing that the tensor can be projected onto the ion density term. The resulting expressions for the force-force time correlation functions all have the mode-coupling theory structure and need to be solved self-consistently. Results of such calculations show that while the hydrodynamic term itself makes a negligible contribution to the friction, the cross-correlation between the electrolyte force term and the HIs term reduces the friction on the ion noticeably, bringing the theory to a much better agreement with known experimental results for several electrolytes. In order to gain understanding of this rather surprising result, we also made an analysis of the length and time scales of relaxation of the various terms. We found that the reason for the relative importance of the cross-term is the overlap of both the length and the time scales of the two processes. Both the ionic interaction and the HIs are of long range. Both relax essentially with the time scale of the Debye relaxation time. In order to further verify the theoretical prediction of the importance of HIs, we carried out Brownian dynamics simulations with and without HIs. Not only do the results agree qualitatively in showing higher diffusion in the presence of HIs, even the agreement between theory and simulation is correct. Therefore, we believe that our prediction of the significant role of HIs via coupling with ion atmosphere relaxation stands confirmed. It is interesting to ask the range of validity of the present scheme. We believe that for 1:1 electrolyte, the present theory with MCT descriptions for both the ion atmosphere relaxation and the hydrodynamic interaction is capable of providing accurate results to 2 M solutions or so. Beyond that concentration, such factors as molecular aspects of ion-solvent interaction and three-body terms can become sufficiently important to deserve special attention and a fully molecular level description, with solvent included, is required. Acknowledgment. We are grateful to O. Bernard for several useful discussions. Note Added after ASAP Publication. This paper was published ASAP on July 8, 2008. The reference section was replaced with a corrected version. The updated paper was reposted on July 25, 2008. References and Notes (1) Wolynes, P. G. Annu. ReV. Phys. Chem. 1980, 31, 345. (2) Onsager, L.; Fuoss, R. M. J. Phys. Chem. 1932, 36, 2689. (3) Barthel, J. M. G.; Krienke, H.; Kunz, W. Physical Chemistry of Electrolyte Solutions; Springer: Darmstadt, Germany, 1998. (4) Fuoss, R. M.; Accascina, F. Electrolytic Conductance; Interscience: London, 1959. (5) Re´sibois, P. M. V. Electrolyte Theory; Harper & Row: New York, 1968. (6) Debye, P.; Hückel, E. Phys. Z. 1923, 24, 305. (7) Chandra, A.; Bagchi, B. In AdVances in Chemical Physics; Prigogine, I., Rice, S. A., Eds.; John Wiley & Sons: 1991; Vol. 80, p 1. (8) Altenberger, A. R. J. Phys. A: Math. Gen. 1981, 14, 957. (9) Ebeling, W.; Rose, J. J. Solution Chem. 1981, 10, 599.
J. Phys. Chem. B, Vol. 112, No. 33, 2008 10271 (10) Friedman, H. L.; Raineri, F. O.; Wood, M. D. Chem. Scr. 1988, 29A, 49. (11) Zhong, E. C.; Friedman, H. L. J. Phys. Chem. 1988, 92, 1685. (12) Bernard, O.; Kunz, W.; Turq, P.; Blum, L. J. Phys. Chem. 1992, 96, 3833. (13) Dufreˆche, J.-F.; Bernard, O.; Turq, P. J. Chem. Phys. 2002, 116, 2085. (14) Dufreˆche, J.-F.; Bernard, O.; Durand-Vidal, S.; Turq, P. J. Phys. Chem. B 2005, 109, 9873. (15) Yamaguchi, T.; Matsuoka, T.; Koda, S. J. Chem. Phys. 2007, 127, 234501. (16) Blum, L. Primitive Electrolytes in the Mean Spherical Approximation. In Theoretical Chemistry: AdVances and PerspectiVes; Eyring, H., Henderson, D., Eds.; Academic Press: New York, 1980; Vol. 5, pp 1-66. (17) Blum, L.; Høye, J. S. J. Phys. Chem. 1977, 81, 1311. (18) Chandra, A.; Bagchi, B. J. Chem. Phys. 1999, 110, 10024. (19) Chandra, A.; Bagchi, B. J. Chem. Phys. 2000, 112, 1876. (20) Chandra, A.; Bagchi, B. J. Phys. Chem. 2000, 104, 9067. (21) Dufrêche, J.-F.; Bernard, O.; Turq, P.; Mukherjee, A.; Bagchi, B. Phys. ReV. Lett. 2002, 88, 095902. (22) Lyklema, J. Fundamentals of Interface and Colloid Science; Academic Press: London, 1995. (23) Jardat, M.; Bernard, O.; Turq, P.; Kneller, G. R. J. Chem. Phys. 1999, 110, 7993. (24) Aldrin Denny, R.; Reichman, D. R. Phys. ReV. E 2001, 63, 065101. (25) Aldrin Denny, R.; Reichman, D. R. J. Chem. Phys. 2002, 116, 1979. (26) McMillan, W. G.; Mayer, J. E. J. Chem. Phys. 1945, 13, 276. (27) Friedman, H. L. A Course in Statistical Mechanics; Prentice-Hall: Englewood Cliffs, NJ, 1985. (28) Lyubartsev, A. P.; Laaksonen, A. Phys. ReV. E 1995, 52, 3730. (29) Lyubartsev, A. P.; Laaksonen, A. Phys. ReV. E 1997, 55, 5689. (30) Ebeling, W.; Scherwinsky, K. Z. Phys. Chem. (Leipzig) 1983, 264, 1. (31) Hansen, J.-P.; McDonald, I. R. Theory of Simple Liquids; Academic Press: London, 1986. (32) Onsager, L. Phys. Z. 1927, 27, 388. (33) Altenberger, A. R.; Friedman, H. L. J. Chem. Phys. 1983, 78, 4162. (34) Onsager, L. Ann. N.Y. Acad. Sci. 1945, 46, 241. (35) Wood, M. D.; Friedman, H. L. Zeit. Phys. Chem. N. F. 1987, 155, 121. (36) Bagchi, B.; Bhattacharyya, S. Mode Coupling Theory Approach to the Liquid-State Dynamics. In AdVances in Chemical Physics; Prigogine, I., Rice, S. A., Eds.; John Wiley & Sons: 2001; Vol. 116, pp 67-221. (37) Sjo¨gren, L.; Sjölander, A. J. Phys. C: Solid State Phys. 1979, 12, 4369. (38) Bhattacharyya, S.; Bagchi, B. Phys. ReV. E 2000, 61, 3850. (39) Bhattacharyya, S.; Bagchi, B. J. Chem. Phys. 1997, 106, 1757. (40) Bagchi, B. J. Chem. Phys. 1991, 95, 467. (41) Bagchi, B. J. Chem. Phys. 1998, 109, 3989. (42) Sung, W.; Friedman, H. L. J. Chem. Phys. 1983, 80, 2735. (43) Turq, P.; Barthel, J.; Chemla, M. Transport, Relaxation and Kinetic Processes in Electrolyte Solutions; Springer-Verlag: Berlin, 1992. (44) Gru¨n, F.; Jardat, M.; Turq, P.; Amatore, C. J. Chem. Phys. 2004, 120, 9648. (45) Ebeling, W.; Feistel, R.; Kelbg, G.; Sändig, R. J. Non-Equilib. Thermodyn. 1978, 3, 11. (46) Kubo, R.; Toda, M.; Hashitsume, N. Statistical Physics II; Springer: Berlin, Heidelberg, Germany, New York, 1991. (47) Kubo, R. ReV. Prog. Phys. 1966, 29, 255. (48) Dufrêche, J.-F.; Jardat, M.; Olynyk, T.; Bernard, O.; Turq, P. J. Chem. Phys. 2002, 117, 3804. (49) Batchelor, G. K. J. Fluid Mech. 1972, 52, 245. (50) Batchelor, G. K. J. Fluid Mech. 1976, 74, 1. (51) Dhont, J. K. G. An Introduction to Dynamics of Colloids; Elsevier: Amsterdam, 1996. (52) Chandra, A.; Wei, D.; Patey, G. N. J. Chem. Phys. 1993, 99, 2083. (53) Ermak, D. L. J. Chem. Phys. 1975, 62, 4189. (54) Jardat, M.; Bernard, O.; Treiner, C.; Turq, P. J. Phys. Chem. B 1999, 103, 8462. (55) Jardat, M.; Durand-Vidal, S.; Da Mota, N.; Turq, P. J. Chem. Phys. 2004, 120, 6268. (56) Mills, R.; Lobo, V. M. M. Self-diffusion in Electrolyte Solutions; Physical Sciences Data 36; Elsevier: Amsterdam, 1989. (57) Kunz, W.; Turq, P.; Bellisent-Funel, M.-C.; Calmettes, P. J. Chem. Phys. 1991, 95, 6902. (58) Cartailler, T.; Turq, P.; Blum, L.; Condamine, N. J. Phys. Chem. 1992, 96, 6766.
JP801796G
