J. Phys. Chem. 1995,99, 10392-10407
10392
Modified Debye-Hiickel Approximation with Effective Charges: An Application of Dressed Ion Theory for Electrolyte Solutions Roland Kjellandert Department of Physical Chemistry, Goteborg University, S-412 96 Goteborg, Sweden Received: January 3, 1995@
The physical basis of the recent dressed ion theory (DIT) by Kjellander and Mitchell (Chem. Phys. Lett. 1992, 200, 76; J. Chem. Phys. 1994, 101, 603) is analyzed and illustrated. In DIT the exact theory for distribution functions in primitive model electrolytes is cast in a form that is virtually the same as in the Debye-Huckel (DH) theory. It is shown in the current work how the utilization of rigorous concepts and exact relationships from DIT (e.g. effective charges of ions and their relationship to the decay length of the pair distributions) can be used to improve the traditional DH approximation and its nonlinear counterpart in a very simple but still substantial manner. The resulting modified versions of the DH theories give values of thermodynamical quantities, effective charge, and decay length (different from the Debye length) for 1:1 and 2:2 aqueous electrolytes that compare well with results from more sophisticated theories and simulations. This constitutes an example-the simplest possible one-of how DIT can be used to improve existing approximations in a new manner.
of this gap. It is a formally exact theory for primitive model electrolytes (i.e. hard sphere ions in a dielectric continuum The Debye-Huckel (DH) theory for electrolyte solutions' is solvent), and it gives a novel insight into the properties of the simplest theory that contains the most prominent feature of electrolytes. Here, "exact" means that the statistical mechanical electrolytes, the screening of charges. It is well-known that theory for the primitive model is formulated without making this theory gives correct limiting laws for electrolytes at infinite any approximations; that is, all physical effects from the dilution and that it is also useful at finite but low concentrations, complete Hamiltonian of the system are included. at least for 1:l aqueous electrolyte solutions. The DH theory In DIT it is shown that the exact theory can be cast in a form is based on the linearized Poisson-Boltzmann (PB) equation. that is virtually the same as in the DH theory provided that The (nonlinear) PB theory has a somewhat larger range of each ion is replaced by a "quasi-particle" consisting of an ion usefulness than the DH theory. A theory that is intermediate with a diffuse charge, a dressed ion. The dress of each ion is between the DH and PB theories is the nonlinear Debye-Huckel an extra charge associated with the ion, and it surrounds the where the average electrostatic potential from the linear central point charge. It is shown in refs 10 and 11 that the DH theory is inserted in the Boltzmann factor of the radial mean electrostatic potential (defined as in the PB theory) distribution function in the nonlinear theory. Following Card satisfies a linear Poisson-Boltzmann equation for the dressed and Valleau," we will call this approach the DHX approximation. ions and that many other aspects of the exact theory are closely Improvements of the DH and PB theories5 range from being analogous to the DH theory. more or less empirical to quite fundamental. Among the latter DIT is not intended as a theory for calculations of electrolyte kinds of improvements one may, for example, mention the properties from scratch, but it can, for example, be used to modified Poisson-Boltzmann (MPB) theory of O u t h ~ a i t e , ~ . ~ analyze data calculated by various approximate theories or by which is quite successful. For some properties, even very simple simulations. Thereby one can extract and reveal several appro ache^^,^ give substantially better results than the DH theory. physically relevant features of electrolyte systems. For instance, Another approach to theories for electrolytes is integral DIT gives a rigorous foundation to concepts such as effective equations, of which the mean spherical approximation (MSA), point charge of an ion and effective permittivity of the solution. the hypemetted chain (HNC) approximation, and their more It also gives explicit formulas for the calculation of these refined counterparts, the generalized MSA (GMSA) and the quantities from the dressed ion charge. reference HNC (RHNC) approximations, are quite commonly The effective permittivity, E , satisfies the condition that the used.2 The HNC and, in particular, the RHNC approximations free energy of interaction between two infinitesimally small and are accurate under a rather wide range of conditions. pointlike test charges dq and dq' separated by the distance r is The accurate theories can only be solved numerically, which given by dq dq' exp(-~r)/[4nEr] when r is large." Here, K is means that much of the physical contents of the results tend to a decay parameter and K - I is a decay length analogous to the be obscured and that the mechanisms in action are hard to Debye length K D - ' . evaluate. The analytically solvable theories like DH and MSA The effective point charge, q*, is the charge with which an are much more transparent in this regard, but they are quite ion is seen to interact with a small test charge dq (as defined inaccurate since they neglect some important effects, and they above) at large separations; that is, the free energy of interaction are therefore unreliable. Thus, there is an annoying gap between between the ion and the test charge equals q* dq exp(-Kr)l accuracy and physical transparency. Recently, a new theory-the [4nEr] when r is large." One finds that q* in general may be dressed ion theory (DIT)'o.''-was formulated, which fills some smaller or larger than the real (bare) ionic charge q, depending on the state of the system.",'* Note that in the DH theory for + E-mail address:
[email protected]. Abstract published in Advance ACS Abstracts, June 1, 1995. point ions q* = q, K = K D , and E = EEO, where EO is the
1. Introduction
@
0022-3654/95/2099- 10392$09.00/0
0 1995 American Chemical Society
Modified Debye-Hiickel Approximation permittivity of vacuum and E is the dielectric constant of the solvent. (As a matter of fact, in DIT this case is used as the reference case in the definition of effective permittivity and charge, so q* = q and E = &EO by definition here.) The effective point charge is related to the decay length K-I by the formulalOsll i
where ni, qi, and &are the density, charge, and effective charge, respectively, of ion species i, k~ is Boltzmann’s constant, and T is the absolute temperature. This formula is an amazingly simple generalization of the definition of the Debye parameter KD given by I
Since qp# qi, it follows that K # KD; in fact, K can be larger or smaller than KD. When the concentration goes to 0, qi+ qi and the standard DH results are recovered. A DIT analysis of 1:l and 2:2 electrolytes in the HNC approximation has been done by Ennis, Kjellander, and Mitchell.’* They made accurate calculations of the variation of K , E, and q* with concentration. Several kinds of transitions between different behaviors of the pair distribution functions were also detected using DIT. Most of these transitions occur at rather high concentrations and will not be considered in the current work. Some of these aspects of the decays of the pair functions have also been treated by AttardI3 and by Leote de Carvalho and Evans.I4 The purpose of the current paper is 2-fold; we will introduce, describe, and analyze dressed ion theory using a minimum of statistical mechanical theory, and-guided by the insights provided by DIT-we will give very simple modifications of the standard DH theory that lead to substantial improvements. This should be of considerable interest in, for example, applied fields such as surface and colloid science, electrochemistry, and biophysics, where simple theories of electrolytes ought to be useful. The main emphasis is, however, to give an example of how the concepts in DIT can be utilized in practice and thereby to outline a new kind of approach that can be further refined in the future. In the usual DH theory for ions with diameter a the ion size is only considered in the interactions with one ion-the central ion-which is, say, of species i. All other ions are treated as pointlike in their interactions among themselves. As shown in ref 10, in the DH theory one can identify an effective charge q i t qr for the central ion,
+
qT= qieKDa/( 1
K,U)
(in DH theory)
(See also section 4.1.) It is only for point ions (a = 0) in the DH theory that qf = qi. Since the surrounding electrolyte solution determines the decay length and since all ions there are treated as point ions and hence have effective charges q,*= qj, it follows that the decay length equals the Debye length KD-’, as predicted by the DH theory. Thus, it is this asymmetric treatment of the central ion and the ions in the ion cloud that makes K equal to KD in the DH theory for ions with size, while, as we have seen, we should normally have K f KD. In the modification of the DH theory presented in this paper, we keep the same identification of effective charge qifor the central ion as in the DH theory, but we maintain the requirement of the exact DIT that the decay length should satisfy K~ = &qiqf/(kBTEEO), where dl ions are treated on the same footing. The self-consistent application of this requirement then leads to the conclusion that all ions have effective charges a** qj
J. Phys. Chem., Vol. 99,No. 25, 1995 10393 and hence that the decay length K - I is different from the Debye length KD-I. The resulting K values for 1:1 electrolytes are equally good as those obtained in the MSA and the linearized MPB theories, and we also obtain quite accurate values of the intemal energy and activity coefficient of the electrolyte solution. To obtain good values of the osmotic coefficient for 1:l electrolytes and to treat the 2:2 case, we make a similar modification of the DHX theory using effective charges. We obtain surprisingly good results for the quite difficult 2:2 system. These modifications of the DH and DHX theories not only give very simple ways to improve the DH approximation but they also give explicit and transparent examples of the DIT analysis of the properties of electrolytes. The modified theories contain some major mechanisms in operation in electrolyte solutions. In section 2 we summarize some basic facts for electrolytes that we need in further developments. Section 3 contains a presentation and analysis of dressed ion theory and its physical meaning and implications. The emphasis is on giving a nontechnical description based on the physical ideas of DIT. The linear version of our modified DH theory is derived in section 4.1, and K and q* are calculated for 1:1 electrolytes. In section 4.2 the nonlinear theory-the modified DHX theory-is derived and calculations are presented for 1:l and 2:2 electrolytes. Explicit formulas for the diffuse charges of the dressed ions are given in the approximations. Calculations of thermodynamical quantities for 1:l and 2:2 electrolytes are presented in section 5, and generalizations of the explicit DH formulas for such quantities are given. (Most of sections 4.1 and 5 can be read before section 3.) Finally, in section 6 some of the results are summarized. 2. Basic Facts
We shall treat electrolyte solutions in the primitive model; that is, the ions are charged hard spheres and the solvent is modeled as a dielectric continuum characterized by a dielectric constant E . For simplicity we shall assume that all ions have the same diameter a, although this is not necessary for our general arguments (see refs 10 and 11). The charge density in the ion cloud around an ion of species j in an electrolyte solution equals
where qi is the ionic charge of species i , ni is the bulk number density, g&) is the radial distribution function of species i and j , r = Irl, r = (x,y,z), and hij(r) = gij(r) - 1 is the total correlation function. In the last equality of eq 1 we have used the electroneutrality condition Xiqini = 0. The origin is placed at the center of thej ion. The total charge distribution, including the charge of this ion, equals
+
+
= qj d‘3’(r) ej(r)= qj d‘3’(r) Cqinjh,(r) (2) i
where ~ 3 ( ~ )is( rDirac’s ) 6 function. The electrostatic potential due to @j’”(r)is called the average electrostatic potential: (3)
This potential satisfies Poisson’s equation: -EEoV 2qj(r)= e?(?-)
(4) where is the Laplace operator. Equation 4 together with the appropriate boundary conditions is equivalent to eq 3 (we demand that yjj(r) 0 when r -). The average electrostatic potential, thus defined, is a central quantity in the Poisson-Boltzmann (PB) equation, which can
-
-
10394 J. Phys. Chem., Vol. 99, No. 25, 1995
Kjellander
be obtained by inserting eqs 1 and 2 into eq 4 and approximating,
g&r) x exp[-Pqiqj(r)],for r 2 a (PB) (5) where /3 = ( k ~ 0 - Iand where (PB) indicates that the expression is only valid in the PB approximation. For r < a we have gV(r) = 0. The exact version of eq 5 is g&r) = exp[-PwU(r)l (6) where w ~ ( ris) the potential of average (mean) force; that is, -Vw,(r) is, on average, the total force on an i ion from the central j ion and its ion cloud. Here, the i ion is one of the ions in the ion cloud, and we take the average of the force for all instances when the ion is located at distance r from the central ion, or, equivalently, we fix the i ion at distance r from the central ion and take the average when all other ions in the cloud are free to move. It follows that in the PB approximation we have assumed that the average force on the i ion equals qiEj(r), where Ej(r) = -Vvj(r) is the electric field originating from the central j ion and its unperturbed ion cloud, which has the charge distribution e ( r ) . Let Fj(r; q) denote the force that acts on a pointlike test charge q placed at distance r from the centralj ion. One can “measure” the electric field Ej(r) by using the relationship
-
F,(r;q) Ej(r) = lim -
q-oq
(7)
The limit q 0 and the vanishing size of the test particle ensure that it does not perturb its environment when the field is measured. Another way of expressing eq 7 is that the force on a pointlike test charge with an infinitesimally small charge dq equals dqE,. If the test particle had a finite charge q or a finite size, Fj would not equal qEj since F, would also contain contributions from the perturbation that the particle induces in its environment; the charge distribution would no longer equal Qj(r). An ion of species i has both a finite charge qi and a finite diameter a. Therefore, the average force on it differs from qiEj, and we see that the assumptions in the PB approximation are not fulfilled; they would only be fulfilled for a point ion with infinitesimally small charge. We ignore in eq 5 that a real i ion affects the ion distribution in its environment, creating its own ion cloud. In other words, the ion-ion correlations in the ion cloud of the central j ion are neglected, a well-known fact. To summarize, in the PB approximation the finite sizes and charges of the ions are fully considered only in the interaction between the central ion and each of the ions in the cloud. The ions in the cloud are treated as pointlike and “cloudless” in their interactions amongst themselves. When calculating the potential of average force acting on an ion of species i in the cloud, it is assumed in the PB theory that the ion behaves as if its charge is infinitesimally small. However, the value of full charge is used in the energy of interaction with the average potential, q@j(r), occurring in eq 5. In the Debye-Huckel (linearized Poisson-Boltzmann) theory it is assumed that each ion in the cloud behaves as if it had an infinitesimally small charge in yet another respect. When qi 0 we have exp(-/3qivj) 1 - /3q$j, so from eq 5
-
-
hij(r)% -Pqiqj(r), for r L a (8) In the DH theory this approximation is used for all values of qi (not only infinitesimally small ones) and eq 8 then replaces eq 5 . Equation 8 is, of course, always a good approximation to eq 5 at large r, where q,(r) is small, but eq 8 is used for all r
2 a in the DH theory. All nonlinear contributions from eq 5 are then neglected. Let us return to the general case. The charge distribution e ( r ) arises from the interactions of the central j ion with all other ions; if these interactions were tumed off, the charge density would be identically 0. (From the point of view of the other ions, such turned off interactions would be the same as if the central ion were not there.) Most, if not all, of the charge density e ( r ) for r L a can be regarded as a polarization charge induced in the electrolyte solution by the electric field from the charge of the central ion, q+Y3)(r),which acts as a source charge. The excluded volume due to the size of the central ion also affects @j(r), making e ( r ) = 0 for r < a and-for nonsymmetric electrolytes-giving some contributions also for r Ia. In the PB theory the last contribution for r 2 a never appears for equally sized ions and @Jr)for r 2 a is solely a polarization charge. Note that the average potential +j(r) is the total electric potential from the source charge 4 i ~ % ~ )and ( r ) the induced charge Qj(r). The change in ion density of species i when we turn on the interaction with the central j ion equals
AnLu)(r)= n,g.,(r) ‘ ‘J - ni= nihV(r) (9) In the DH theory only a linear response to the total electric field is considered, and we have from eqs 8 and 9 Any)(r) -Pniqiqj(r), r a (DH) (10) where (DH) indicates that this is only valid in the DH approximation. The induced polarization charge then equals i
I
Introducing the Debye parameter KD defined as
we obtain
For point ions ( a = 0) we have
@ F ( r ) q,13‘~)(r)- EE,y$q,(r)
(DH)
(13b)
Note that the response to v J ( r )in eqs 10-13 is both linear and local; the induced charge at one point is proportional to the potential at the same point. Since v J ( r )originates from both the source and the polarization charges, the response is selfconsistent (within the DH theory). Equations 4 and 13 yield the linearized Poisson-Boltzmann equation. For point ions this equation is
-EEOV2v,(r) = ~ , 6 ‘ ~ ’ (-r )&&,&qJ(r) (DH) (14) (Usually it is written without the source term, which vanishes for r > 0.) In what follows we will extend the DH treatment of electrolyte solutions, formulating an exact theory, the dressed ion theory, and we will give essential improvements of the approximate DH theory. An important ingredient is the definition and use of effective quantities, such as effective charges of ions. To define such quantities, we will need a standard case with which to compare. The simplest possible case-the DH theory for pointlike ions, exemplified by eqs 13b and 14-will prove to be suitable as such a standard. Therefore, we will compare with this theory when treating the general case.
Modified Debye-Huckel Approximation
J. Phys. Chem., Vol. 99, No. 25, 1995 10395
3. The Physical Basis of Dressed Ion Theory
In the previous section we have seen that the DH theory gives an approximate but self-consistent linear response of the electrolyte solution due to the interaction with a central ion. We have also seen that the particular formulas used in evaluating this response would be correct if all other ions were pointlike and had infinitesimally small charges. To generalize the DH theory, it is fruitful to consider how a bulk electrolyte solution with real ions responds to an electric field originating from some source charge, which, for example, can be a central ion. Since we will consider an infinitely large system, we may without loss of generality regard this source charge as extemal to the system (Le. the system consisting of all other ions). Initially, we will restrict ourselves to the linear response regime and consider the charge distribution induced in the electrolyte solution by a small electric potential 6Yext(r),where the superscript “ext” denotes that the potential originates solely from the extemal source charge. The change in extemal potential from 0 to dYext(r’) for all r’ causes the density of ion species i at point r to change from n, to nl 6n,(r). It follows from linear response theory (see Appendix A) that
+
- r’l) +
6nl(r) = -BnlJdr’[qld‘3’(lr I
Using eq 2 we can write eq 15 as
e,“’(Ir - r’l) GVxt(r’)
(16)
The charge density of the electrolyte changes from 0 to 6g(r) = Zlqr6nl(r); that is, &(r) is the polarization charge induced by the extemal field. From eq 16 follows
(22)
In fact, as shown by Kjellander and Mitchell,’’ we have (cf. eq 2)
where the function huo(r) can be calculated from hu(r) in a straightforward manner (see Appendix A). Except at high concentrations, huo(r) has shorter range than h&) and &(r) shorter range than pi(r).”x12 It follows from eq 22 that
= ldr’Bo(lr - r’l) dY(r’)
(24)
where
Thus P(r) may be denoted as the “self-consistent, linear response function of p~larization”.’~It is self-consistent since it gives the polarization response due to the potential that partly arises from the very same polarization. Compare eqs 24 and 25 with eqs 17 and 18. By comparing with eqs 19 and 20 we see that the response functions in the DH case can be obtained by approximating
e%) = qr 6‘3’(r) (DH)
(26)
Le. by putting all h$r) equal to 0. Thus, in the DH theory
6g(r) = - B & z , S d r ’ g,fof(lr - r’l) dWxt(r’) I
= Jdr’ B(lr - r’l) d\Y“”‘(r’)
dni(r) = -Bn,Jdr’ e’’(lr - r’l) SY(r’)
I
&n,h,~(lr - r’b] dpxt@’) (15)
6nl(r) = -BnlJdr’
change in density is obtained from the total potential by the relationship
B&)(T) = -&E&; (17)
where
d3’(r)
(27)
while the correct P ( r ) is nonlocal; there are short-ranged contributions from h:(r) present:
Bo(r)= -E&&D
2
6(3)(r) - B&,q,nfn,hI,O(r)
(28)
Y
is the linear response function of polarization. We may compare eqs 16 and 17 with the DH formulas 10 and 11, which we may write
AnY)(r)= -Bn,Jdr’ qi 8(3)(lr- r’l) q j ( f )(DH) (19)
ej(r)= -BCqiniJdr’ qi 6(3)(lr- r’l) qj(r‘)(DH) (20) i
The linear response in the DH case is, as we have seen, entirely local, which in these equations is expressed by the Dirac 6 function. The response functions in eqs 16 and 17 are nonlocal. The most important difference is, however, that the potential v,(r) in eqs 19 and 20 is the total potential due to the source charge and the induced polarization pj(r). The, potential 6Yext(r)in eqs 16 and 17 is solely from the source charge. To generalize the DH response function, we must in the general case find the linear response function that gives the response due to the total electric field. The total potential from the extemal source charge and the induced polarization is
In analogy to eq 16, define the function e&r) such that the
Let us now apply eq 24 to the case where the source charge is an ion located at the origin. Since a real ionic charge is not small, the electric potential is so large that nonlinear contributions to the polarization response in general are important in the neighborhood of the ion; that is, there are contributions in addition to the linear one given in eq 24. Therefore, we shall fvst treat a case where the nonlinear contributions are negligible before we do the general case. We assume that the central ion has a fractional charge, 6% (in principle it should be infinitesimally small), so the source charge distribution is 6% d(3)(r). The electrostatic potential from this source charge gives rise to a polarization charge density, which we denote de,(+ The total electric potential from the source charge and the polarization charge will be denoted dv,(r), and it can be obtained from the charge distributions by applying Coulomb’s law (cf. eq 21 with d Y = dv, and 6g = de,). From eq 24 follows that the total charge density then equals
+ = dq, d(3)(r)+ J dr’ Bo(lr - r’l) dq,(r’)
dg?(r) = 64, 6‘3’(r) de,@)
(29)
Note that eq 29 refers to the situation when only the central ion has a small charge; all other ions have full charges and they interact and correlate fully with each other.
10396 J. Phys. Chem., Vol. 99,No. 25, 1995
Kjellander
We may compare eq 29 with the DH approximation for point ions (Le. our reference case), eq 13b, which can be used for a central ion with either an infinitesimally small charge or a full charge. In both cases only the linear response is included in the DH approximation. Using the DH response function in eq 27, we can write eq 13b in the same form as eq 29:
+
= qj 6(3)(r) Jdr’ B&)(lr - r’l)
vj(r‘>(DH)
(30a)
Altematively, we can write eq 30a as (cf. eq 20)
@y(r) = qj CY3)(r)pxqiniJdr’ qi 6‘3’(lr- r’l) Vj(r‘) (DH) (30b)
a diffuse contribution (the last term), which is a short-ranged charge distribution that surrounds the central point charge of the ion. This charge distribution is a well-defined part of the ion cloud distribution around the ion. The rest of the ion cloud distribution equals B#(r); see eq 34a. As mentioned above, eF(r) has a shorter range than e,(r) (except possibly at high concentrations). The effective diffuse charge of each ion can therefore alternatively be regarded as the charge of a “dressed ion”, a quasi-particle that we may use to replace the bare ion when formulating the theory. To see that the resulting theory for the dressed ions is virtually identical to the DH theory, let us combine eqs 4 and 34 to obtain the exact version of the linearized Poisson-Boltzmann equation:
i
We now proceed to the general case of a central j ion with full charge qj. Since we want to formulate a generalization of the DH theory, it is desirable to keep the structure of eq 30a,b:
+
total charge = source charge polarization charge (linear response) (31) Equation 29 for a fractionally charged central ion also has the structure of eq 31, but we want to have the corresponding formula for a fully charged ion. The total charge ept(r)is given by eq 2. However, @,(r)in eq 2 contains in the general case both linear and nonlinear response contributions, so we have to separate out the linear part, which we denote e&). It is given by (cf. eq 24)
@,”(r)= Jdr’Bo(lr - r’l) V j ( f )
= -PxqiniJdr’ e;(Ir
- r’~)V j ( { )
-eeoV21C;(r)= @;(r)
+ Jdr’
Bo(lr - r’l) V I ( / ) (35a)
The main difference from the ordinary DH equation 14-apart from that eq 35a is exact-is that the linear response term is nonlocal, which is due to the diffuse nature of the dressed ions. If we write eq 35a as
-EEoV2VJ(r)=
- pxqJnJJdr’@:(Ir
- r’l)
$I:(”)
I
(33)) we see that the dressed ion charge acts as both the source of the average potential and the charge that interacts with this potential; the integral in the last term equals the energy of interaction. Let us define
(32)
i
If we identify e,#(r)with the last term in eq 31, the term “source charge” cannot, of course, equal the true source charge, q, d(s(r). In this term we must also include the difference ej(r) - @,#(r)which contains the nonlinear contributions, so eq 31 reads
+
+
= [qj6(3)(r) ej(r)- @;@)I #(r) (33) The sum inside the square bracket thus plays the role of an effective source charge when we insist on including only the linear response in the term “polarization charge”. The amazing fact is that the square bracket exactly equals &r), i.e. the same function as in eqs 22 and 32, as shown by Kjellander and Mitchell.’o,” This means that we have @?(I.)
= ejO(r) - B&n,Jdr’
@;(b- r’l) Vj(r‘) (34a)
i
which has exactly the structure of the DH expression 30b, the only difference is that in the DH case we have replaced ejo(r) with q, d(3)(r),i.e. the approximation stated in eq 26. Using eq 25 we can write eq 34a as
@y(r) = e;(‘) + l d r ’ Bo(lr - r’l) Vj(r‘)
so we have from eq 32 (37) From eqs 8 and 26 we see that in the DH approximation for point ions hu(r) = hif(r). In the general case we
+
hu(r)= h J r ) h,#(r) (38) and again it follows that the DH approximation is obtained when we set huo(r) = 0. We may call h / ( r ) the “DH-like part” of hu(r). As mentioned above, huo(r)has shorter range than h&) except at high concentrations, so eq 38 constitutes a split of hu(r) into a short-range and a (relatively) long-range part. As we have seen, the exact formalism for the average potential and the distribution functions takes the same form as in the DH theory when we use dressed ions instead of bare ones. This does, of course, not mean that the results of the DH theory remain true: the solutions of, for example, eq 35 are, in general, quite different than those of the usual DH equation. However, some aspects of the solutions are similar in the two cases. In the DH approximation for point ions the average potential decays as
(34b)
which should be compared with eq 29 and the DH expression 30a. Note that eqs 34a,b are exact: no approximations have been made. Thus ey(r) plays the same role in the general case as the bare ionic charge plays in the DH approximation. It can be interpreted as an effective diffuse charge of each ion, a renormalized charge. If we use this effective diffuse charge instead of the bare charge qj d(3)(r)of each ion, the exact theory accordingly assumes the same form as the DH theory. From eq 23 we see that e?(‘) consists of the bare point charge plus
(39) where KD is given by eq 12. In the exact, general case the average potential decays for large r aslo,ll
where q,* is an effective point charge that, in general, is different from qj, and E is an effective permittivity that generally is different from the permittivity EEO of the solvent. The decay
J. Phys. Chem., Vol. 99, No. 25, 1995 10397
Modified Debye-Huckel Approximation parameter K satisfies the relationship
which is a simple generalization of the definition of K D in eq 12. We see that K K D in general. Distinguish the effective point charge %* from the dressed ion charge (effective diffuse charge) &r). They play different roles in the theory but are closely related; in effect, qj* is the charge that a point ion must have in order to interact like a dressed ion, which has a somewhat spread out charge. More precisely, when the effective point charge and the dressed ion are placed at the same position in the large r tail of the average potential, they have the same interaction energy," Le.
*
--
-
when r S d r ' @:(lr - r'l) Vj(r') qi*Vj(r) Inserting eq 40 in this equation, one can show that"
qi* = e J w d r e:(r)r K
(42)
for large r. The interaction free energy between that charge and another small point charge dq' separated by a distance r thus decays as
-
dq dq' e-Kr
(45) 4nEr when r is increased. Note that the actual charge dq occurs in eq 44, while the effective point charge %* appears in eq 40 for an ion with a full charge. This difference can be traced back to the Occurrence of the infinitesimal charge as the source charge in eq 29, while the dressed ion charge is the source in eq 34b. The formula corresponding to eq (45) for two fully charged ions is'0,'I wij(r)
-
qi*qj*e-Kr
4nEr
(46)
for large r, where wu is the potential of average force (eq 6). Thus two ions in an electrolyte solution interact as if they had another charge than their real (bare) charge. At large separations they interact with apparent (effective) point charges given by g*, I = ij. The value of q1* depends on the state of the electrolyte solution. When the concentration goes to 0, q1* 41, K K D , and E &EO, so the ordinary DH theory results are recovered. Since eq 6 yields gu(r) 1 - fiwu(r) for large r, where wu is small, it follows from eq 46 that
-
-
where K = Re(K), o = I ~ ( K and ) , Qj*and q5j are real parameters (Re denotes the real part and Im the imaginary part). The wavelength 2do and the decay length 1/Kare common for all species, while Q,* and the phase shift 4, are species dependent. For hu(r) one obtains
sinh (Kr)
which, in fact, can be taken as the basic of qi*. Incidentally, eq 43 can be written qi* = QF(iK), Le. the Fourier transform of &r) evaluated at iK (where i is the imaginary unit). The effective permittivity E and the decay parameter K are properties of the electrolyte solution. For example, the average electrostatic potential dV(r) due to an infinitesimally small point charge dq placed in the solution is"
dw(r)
have seen, the function huo(r) in eq 38 has a shorter range than h&) (except at high concentrations). Therefore, h$'(r) also decays according to the right-hand side of eq 47. This also follows from eq 36 combined with eqs 42 and 40. Equations 40 and 44-48 actually have a limited range of validity; they are not valid at high concentrations. At a certain concentration, which depends on the system, the parameter K turns complex (and so do q1* and E ) . Physically, this means that the decay becomes damped oscillatory rather than monotonic for large r,
-
-
where the phase shift 8u is real. This kind of change in behavior was first found theoretically for electrolyte solutions by KirkwoodI6 and has since been found in most theories of electrolytes. Other kinds of changes in behavior for w&r) and h&) also occur at high concentrations; see refs 12-14. This happens when the decay length of huo(r) becomes longer than l/Re(K), and eq 47 or 50 therefore do not give the dominating contribution to hu(r) at large r. Here we shall not consider these complications any further. Before closing this section we shall make some connection to the nonlinear Poisson-Boltzmann theory. From eq 6 follows that -@wij(r) = ln[l hu(r)] for r 2 a. Defining &(r) = hi)'(r) ln[l hij(r)] - hij(r) and using eq 38, we have
+
+
+
-Pwij(r) = h,#(r) +x,(r)
for r
a
(51)
At large r, where hu(r) is small, we can expand the logarithm infi, in a MacLaurin series and obtain Jj(r) x h,O(r) - [h,(r)12/2 (52) Except at high concentrations, we see thatJj(r.0 decays quicker to 0 than hf(r) when r = (in accordance with what we already know). Furthermore, when the concentration goes to 0, we have" hvo(r) % [hu(r)I2/2which meansfl(r) becomes even smaller. In fact,fij(r) 0 when the concentration goes to 0, while hu#(r) stays finite (see Appendix A), so -@wij(r) % hf(r). Inserting this approximation in eq 6 and using eq 36, we obtain
-
-
g&r)
= Io,exp[-Pjdr'
@:(lr-rfl)
~,(r')],r L a
(53) r KD when a > 0. This result may appear quite striking, considering that the usual DH theory predicts that K = KD but otherwise has exactly the same formulas for e, and lyj. However, in the DH theory the central ion is treated in a way different from the rest of the ions. Except in the interaction with the central ion, the ions in the ion cloud are treated as if they behave like point ions with infinitesimally small charges, and, as we have seen, for such ions qi* = qi. Since K is a property of the electrolyte solution surrounding the central ion, it follows that K = KD in the DH case. It is only the central ion that has an effective charge q,* t q, in the DH theory. In contrast to this, in our treatment here we have made no assumptions that the central ion is diTerentfrom the other ions and, as a consequence, K f KD. All ions of species 1 have effective charge ql* f 41, 1 = 1, 2, .... It follows that the usual DH theory has a kind of built-in asymmetry in the pair distribution functions, which may not appear obvious considering that from eqs 8 and 58 (with K = K D ) follows that hij(r) = -
0 when
These results are identical to those of the usual DH theory for ions of diameter a, except that we have not yet specified K (in the DH case K always assumes the value K D ) . If, instead of the central j ion, we have a pointlike, infinitesimally small charge Qq placed at the origin in the electrolyte, the resulting average potential &#(r) can be obtained from eq 58 with a = 0:
pqiqIeKDae-KDr
r l a(DH)
+
~x&E 1 ~ K( ~ u ) ~ ‘
(61)
an expression that is symmetric in the indices i and j . However, using the identification (60) (with K = KD), eq 61 should be written pq .q.*ePKDr hij(r) = - 4 n ~ & , r , r 2 a (DH) (62) and we see that the symmetry of eq 61 is rather accidental. In contrast to this-provided the pair distribution has the functional form hij(r) = aiajF(r) (where aiand aj are constants and F(r) is some function)-eqs 41, 57, and 60 imply that h, must in our present approach equal (see Appendix B for proof)
’
pqi*qj*e-Kr ,rLa 4ns&,r which is symmetric in q*. This is an important difference from the ordinary DH theory. Thus, eq 61 is replaced by hij(r) = -
(59) If we compare with the general formula (44) we conclude that E = EEO for the electrolyte in the approximation given by eq 55. Comparing eqs 57 and 58 with the general results 48 and 40, taking E = EEO into account, we can identify the effective point charge of the j ion as qjeKa qj* = 1 KU and we see that qj* > q if a > 0. Note that the only effect we have included is the void of radius a that the j ion creates in
+
pqiqje2Kae-Kr ,r2a (64) 4nc~,(1 Ka)*r in the current approach. (Note that eq 63 is a direct consequence from eqs 69 and 70 below, so the assumption about the functional form is not extraneous in the approach.) If we insert eq 60 in eq 41 we obtain an equation for K , h&r) = -
+
+
K ~ / K= , ~eKa/(l
(654 which can be solved numerically. Introducing the notation Y KU)
Modified Debye-Huckel Approximation
9 .x
4
-
J. Phys. Chem., Vol. 99,No. 25, 1995 10399
..
4
m
..
cv,
. Y
'0.8
1.0
1.2
1.4
.
Y v
1.0
0.5
0.0
1.5
D' a Figure 1. The section of each curve below the symbol shows the decay
parameter K (inverse decay length) relative to the inverse Debye length plotted as (K/KD)* as a function of K D for ~ a 1:1 electrolyte solution from various theories. (Note that KN is proportional to the square root of the concentration.) The plot also shows the ratio of the effective point charge q* of the ions to the ionic bare charge q, since q*/q = ( K / K D ) ~for symmetric electrolytes. The symbol indicates the point where a second solution of the equation for K and q* (section of each curve above the symbol) merges with the first, dominating solution. At higher concentrations the solution becomes complex and the decay of the pair distribution oscillatory (see text). The curves show the results from MSA (- * -), HNC and our modified DH approximation (-). The HNC data are taken from ref 12 and show results for a system of monovalent ions with a = 4.2 A, E = 78.358, and T = 298.15 K. The inset shows sections of the curves close to the symbols on a different scale. KD,
(.*a),
= KIKDand z = KM, it can be written
+
v2 = ev7/(i vt)
(63)) and we see that KIKD is a function of KDU alone. Equation 65b has two real, positive solutions, one of which goes to 1 when z 0; that is K K D when the concentration goes to 0 (or when a 0). This is the wanted solution. The other solution goes to infinity when z 0. At z = 1.35 the solutions tum complex; at this point the second solution merges with the first one. In Figure 1 the two solutions are plotted and compared with the corresponding plots from the mean spherical approximation (MSA) and the hypemetted chain (HNC) approximation for a 1:l electrolyte in water (the latter taken from ref 12). The linearized version of the modified PoissonBoltzmann theory (LMPB) by Outhwaite6*20gives results for K ~ I Kthat D~ are very close to those of MSA and are not shown. In all cases K values turn complex when the second root merges with the fiist; this point is marked with a symbol in Figure 1. For symmetric electrolytes we have from eq 41 K ~ / K=Dq*/q, ~ so Figure 1 also shows how the effective point charges of the ions vary with KN in this case. In MSA v = K / K D is the solution2' of vt(vt - 2p) = 2p2(eYr - l), where p = [(l 2 ~ ) " ~1112, and we see that KIKDis a function of KDU alone also in this case. The latter is true also in the LMPB theory. The HNC approximation gives more accurate results. In this approximation, as well as in the exact theory, KIKD is not only a function of KW but, for symmetric electrolytes, it depends also on #?*= # ? q 2 / ( 4 m ~(or ~ )any two combinations of these dimensionlessparameters). When dealing with nonsymmetric electrolytes, there are even more independent parameters specifying the system. The HNC case in Figure 1 has #?*= 1.704, which corresponds, for example, to a 1 :1 electrolyte with a = 4.2 A when E = 78.358 and the temperature is 25 OC. The dependence of KIKDon /?* implies that the curves showing KIKD versus KW
- -
-
-
+
cannot coincide for electrolytes with different a values (at constant q, T, and E). For the HNC case in Figure 1 this dependence is small; on the scale of the figure the corresponding curve for a = 5 8, is hardly distinguishable from that for 4.2 A; see Figure 3 and ref 12. Thus, the #?*independence of the results from our current theory and from MSA does not matter much in this regime. Considering how simple OUT approach is compared to the MSA, LMPB, and HNC theories, it does remarkably well. In fact, it gives at least equally good K values as the MSA and LMPB theories. However, it is not realistic around the transition to complex K values and beyond, since in eq 55 we have explicitly assumed that there is only one value of K . (To allow for more than one K , one must include more exponential terms.) Nevertheless, we can conclude that the main mechanism that gives rise to the increase in KIKD and q*/q with concentration is included in our approach. This mechanism is the simple exclusion of ion cloud charge in a radius a around the ions, which leads to a buildup of appreciable charge density at larger r values than for a point ion. It is essential that this exclusion is considered for all ions. If we would only have an exclusion region around the central ion, we would get the usual DH result and KIKD= 1. Note that the qualitative behavior shown in Figure 1 was already obtained by Kirkwood16 when he included effects of exclusion around the central ion as well as other ions. His equation for v = K / K D is v2 = cosh(vt), and K t u m s complex at KDU = 1.03. The exclusion mechanism that leads to K f KD is included in all theories that consider the size of all ions in the electrolyte. A virtue of our approach is that it demonstrates the effect of exclusion in a particularly transparent manner. Since the linear theories (Le. all mentioned above except the HNC theory) do not include any dependence on parameters other than KDU, their behavior must be qualitatively wrong in some regimes. For instance, they all incorrectly predict that KIKD 1 when a 0 and that we always have K/KD > 1 independently of /?*. In fact, when K D becomes sufficiently small, KIKDtums virtually independent of a, and for symmetric electrolytes KIKD eventually becomes less than 1. This follows from the exact, asymptotic behavior of K ~ / K Dfor~ small values of the coupling parameter A = K D ~ B(where ZB = ,8e2/{4m~o}is the B j e r " length) derived by Kjellander and Mitchell:"
-
-
-K2 =I+2
KD
A In 3 4
(z,
+ ZJ2 + In A (z12 6
+ zlz2 + z):~ + O(A2) (66)
where zi = qile is the ionic valency. The ion size dependence enters first in the higher order terms (A2 and higher); the terms shown here are of entirely electrostatic origin. For symmetric electrolytes, where zl = -22, the second term vanishes and the A2 In A term, which is negative when A < 1, dominates at small A. (Note that A = KDU /3*/z2 in this case.) For the 1:l electrolytes in Figure 1 these effects have a small influence, but they can nevertheless be seen in the low KDQ regime shown in Figure 2, where the HNC curve clearly goes below 1 at low KDU. (The HNC theory satisfies eq 66.) On the other hand, for 2:2 electrolytes in water these effects are much more prominent and there is an appreciable ion size dependence of the KIKDversus KDU curves.12 Then the linear theories give poor results. The same is true for 1 :1 electrolytes with small a values, since an increase in z is equivalent to a decrease in a (B* increases in both cases) provided the concentration is adjusted so KDU is the same. In these cases it is essential to use nonlinear theories. Our present approach with effective charges can readily be carried over into a nonlinear
Kjellander
10400 J. Phys. Chem., Vol. 99, No. 25, 1995
1.30
I
’
I
I
I
1
1.25
L
I / .
c
.
Y
i.e. inside a spherical hole in a radially symmetric charge distribution with a point charge qi at the origin, equals
’
1.10
1.00
1
From eqs 37, 69, and 73 follows that @j#(r),the Debye-Hiickellike charge distribution, is nonzero inside the core. This is due to the fact that @j#(r)does not explicitly consider the effect of the hard core interaction but only the electrostatic potential. Inserting eq 73 in eq 72 and using eq 58 at r = a, we obtain from eq 23
1
L L L J .....
0.0
0.2
0.4
0.6
0.8
‘D a
Figure 2. Same plot as in Figure 1 shown in the low
KM
regime.
theory which will be presented in the next section, and we shall wait with the treatment of the 2:2 electrolytes until then. We shall finish the treatment of the linear theory with an investigation of the dressed ion charge distribution, @(r) implied by the current modified DH approximation. The relationship to the effective point charge q* will be investigated. We shall also see how the approximation in eq 55 can be motivated from a more general point of view. This gives a suitable starting point for the development of the nonlinear theory later. The basic approximation we use is (cf. eq 51 withfij 0 and eq 36)
-Pw,(r) m h,#(r) = -BJdr’ e:(lr - r’l) qj(r’), r 2 a (67) Le. the same as in eqs 53 and 54. Equation 53 can be written
From eqs 38 follows hIF(r)= gI,(r)- 1 - hI,#(r). Thus, eq 68 implies that we have expressed hLF(r)and hence @?(r)(see eq 23) in terms of hI,#(r). Furthermore, we use the large r behavior of the integral in eq 67 (see eq 42) as an approximation for all r values, Le. h,#(r> -Bq,*q,(r) for all r (69) (Note that this does not imply that we have taken @?(r)= ql* 13(~)(r); we have only approximated the integral in eq 67, while Q?(r) for all r is determined by h,#(r) from the relationships above.) In the linear theory we linearize eq 68 and obtain
h&r)
=
-1, r
The (osmotic) pressure of the solution equals22
where the first term is the ideal pressure Pideal.Due to the discontinuity of the pair potential at the hard sphere diameter, there will be one contribution to the integral from r = a (the contact term) and one from r > a. For Coulombic interactions the latter is proportional to the excess energy. Defining the osmotic coefficient from 4 = P/Pdea',we obtain
where the second term is the contact term. The excess chemical potential p:x and hence the activity coefficient yi of species i can be obtained from a coupling parameter integration,
G
Physically, what happens is that the @p(r) contributions outside the hard core region of an ion originate from the nonlinear part of the response of the electrolyte due to the interaction with this ion. If we want to mimic the total response in a model with only linear response (our reference system), we must change the charge of the ion, i.e. use a q* # q, to get equal values of the response at large r. As regards the dressed ion charge distribution inside the hard core, its physical role is the same as in the linear theory (see previous section). Note that the reason for the initial negative deviation in q* and K is the same as that giving rise to the general, asymptotic
where 5 is the coupling parameter, uij(r;E)for 0 < 5 < 1 is the pair interaction between a partially coupled ion of species i and the other ions (the latter interact fully with each another), and n.,gI1(r;5) .. is the density of species j around the partially coupled i ion. The potential ub is a continuous function of E such that uij(r;O)= 0 and ui,(r;l)= ui,(r) for all r. Note that the coupling involves a change in both size and electric charge of the ion; both the size of the hard core region and the ionic charge are 0 when 6 = 0, and the ion diameter is a and the ionic charge is qi when 6 = 1. The coupling may, for example, be performed
J. Phys. Chem., Vol. 99, No. 25, 1995 10403
Modified Debye-Huckel Approximation in such a way that the size of the uncharged ion is fust increased to its full size before the charge is increased; this is a coupling procedure in two stages. To be precise, let us fust increase the radius of the hard core exclusion region around the uncharged “ion” (i.e. the distance of closest approach of other ion centers) linearly from 0 to a and then increase the ionic charge linearly from 0 to qi. Let n,gij[A*Qi](r) denote the density of j ions around an ion of charge Qi in an electrolyte solution of fully interacting ions. The ion has a hard core exclusion region of radius A. Using the two-stage coupling procedure, we can readily calculate yi from eq 87 and therefore the mean activity coefficient coefficient y+ (defined from n In y* = 2ini In yi). It is straightforward to show that
The corresponding formula in the DH theory is
The second term is the contact term, which often is left out from the DH formula for 4. However, it should be i n ~ l u d e d , ~ as evident, for example, from simulation studies of the primitive modeL4 The mean activity coefficient in the modified DH theory, as obtained from the coupling parameter calculation in eq 88, equals 4na3n -
In
?’* = 3
K ~ * K
8nn( 1
+
KU)
(94)
and the corresponding DH result is 4na3n
3
KD
In y* = -3 8nn(l The first term is a contact term; note that n,gij[A,o](A)is the contact density of j ions at the surface of an uncharged hard sphere immersed in the electrolyte (this density is actually independent of index i since Qi is 0 here). The mean activity coefficient can altematively be obtained from 4 by a thermodynamic integration from zero concentration to n. Since 4 - 1 is proportional to Jn at low concentrations, a suitable standard formula is
which relationship is a consequence of the Gibbs-Duhem equation. There is no guarantee that the values of y* calculated from eq 88 and from eqs 89 and 86 are equal when approximate distribution functions are used. This kind of inconsistency-Le. that the value of a quantity depends on which route one takes when it is calculated-is normal in approximate theories. In fact, it is impossible to avoid such inconsistenciescomplete1yi.e. for all possible quantities-unless exact distribution functions are used. For sophisticated approximate theories there are, however, only small differences, if any, between the values of each relevant quantity calculated by different routes. Let us now evaluate the various thermodynamical quantities in the linear version of our modified Debye-Hiickel theory. Using the distribution function from eq 63 and the relationships 41 and 65a, we obtain from eq 84
8nn(l
+~ a >
(90)
where UeVN is the excess energy per particle. Note that K D 2 h is constant, independent of the concentration. Equation 90 should be compared to the corresponding result in the ordinary DH theory:
Likewise, the osmotic coefficient in the modified DH theory is obtained from eq 86: #=1+--
ka3n 3
K
24nn(l
~
+
K KU)
(92)
+
K,U)
(DH)
(95)
The first term is the contact term discussed in connection with eq 88 above. It is often left out from the DH expression for the activity coefficient, but should clearly be present for the primitive model; it equals the change in free energy when creating a hole of radius a, where the centers of the surrounding ions cannot enter. In both the DH and modified DH approximations gij[A*o](r) for r 2 A and evaluation of the first term in eq 88 gives the first term in eqs 94 and 95. It is common to assume in the DH theory that the electrolyte solution becomes ideal when the charge of all ions are set equal to 0, and then the first term of eq 95 should not be present. However, its inclusion gives a much better performance of the DH theory (see below). When In y+ is calculated in our modified DH theory from the osmotic coefficient via eqs 89 and 92, a result different from that in eq 94 is obtained:
This is an example of the common inconsistency of approximate theories mentioned above. The corresponding DH result from eqs 89 and 93 is
(97) and is different from eq 95 for the same reason. The contact term here originates from the contact term in eq 93 and is identical to that in eq 95. The formulas for thermodynamical quantities from our modified DH theory and the normal DH theory are very similar, eqs 96 and 97 being an exception. We expect that the linear theory is useful for 1:1 electrolytes, while nonlinear effects are clearly important for 2:2 electrolytes. Therefore, we shall apply our present results to 1:1 electrolytes, while the 2:2 case will be treated using the DHX theories later in this section. In Figure 7 we have plotted the UeX/NkBTresults from eqs 90 and 91 for a 1:l aqueous electrolyte, and also included are the corresponding Monte Carlo (MC) simulation results taken from ref 4. Our modified DH theory gives excellent results and constitutes a clear improvement compared to the normal DH theory. The deviation of the latter is, however, not very large. The vertical line in Figure 7 shows where the pair correlations become oscillatory (as obtained in the HNC theory12). At higher concentrations neither the modified DH
10404 J. Phys. Chem., Vol. 99, No. 25, 1995
Kjellander
0.0
-0.1 Em
-0.2 -0.3
I
-0.4 -0.5
-
‘ ‘.
DH lim I
,
/
‘
I
9’
-0.1
F
1.0
MDH /
]
h
L
0.0
0.2
0.4 c1/2
Figure 7. (Reduced) excess intemal energy per ion UeXINkeTas a function of the square root of concentration c (measured in M) for a 1:l electrolyte of ionic diameter a = 4.25 A when E = 78.5 and T = 298.15 K. The curves are from DH theory (. *.), labeled “DH’,modified DH theory (- - -), labeled “MDH’, and MC simulations (B). The DH limiting law is shown as (-), labeled “DH lim”. The vertical line shows
the concentration where the pair distribution function becomes oscillatory according to the HNC theory. The MC simulation data are taken from ref 4. nor the normal DH theory is realistic since both give monotonic decays. As we have seen, below the point where the oscillations start, there is a second K root, K’, present. It gives rise to a further, exponentially decaying term in the correlation function. This additional term gives important contributions at least when K’ w K (Le. at concentrations near the point where oscillations start), and it is absent from the present DH theories. Therefore, one expects deviations from simulation and HNC results even before the oscillations start. Accordingly, the modified DH theory-in its present form with only one exponential term-is unreliable at concentrations higher than those shown in Figure 7, and at increasing concentrations it eventually breaks down. To extend this theory to the oscillatory regime, it is essential to have both exponential terms present: it is their sum that turns into the oscillatory function in eq 50. As shown by Card and Valleau4 the DH theory does not give particularly good values of the osmotic coefficient (from eq 93) for 1 :1 electrolytes, and it is necessary to go to the DHX theory to obtain good values of 4. The same is actually true in our modified DH and DHX theories, so we leave this until the discussion of the nonlinear theory below. The mean activity coefficient for a 1:l electrolyte obtained from eqs 94-97 is shown in Figure 8 together with the corresponding grand canonical Monte Carlo (GCMC) results of Valleau and C ~ h e n .As ~ ~expected from the rather poor values of 4 in the linear theories, the activity coefficient from thermodynamic integration of 4, eqs 96 and 97, is not very good (curves labeled MDH(4) and DH(4)). The results from the coupling parameter integration, eqs 94 and 95, are, however, quite good (curves MDH and DH 2). Note that the inclusion of the contact term in eqs 94 and 95 is essential for this. The performance of the usual DH expression for In y*, i.e. eq 95 without the contact term, is quite poor (curve DH l ) , a wellknown fact. The difference in performance of the DH and the modified DH theories, eqs 94 and 95, is not large, and due to the errors involved in the GCMC data, it is not easy to determine whether the modified DH is better or not, but up to about 0.6 M (cl/* 0.8) the latter theory apparently performs slightly better than the ordinary DH theory. At somewhat higher concentrations the pair distribution becomes oscillatory (cf. discussion above), so the results of the DH theory there are not very relevant. Let us proceed to the nonlinear theory, the modified DHX theory. The excess intemal energy and the osmotic coefficient
0.6
0.8
1.0
/ M1/2
Figure 8. Logarithm of the mean activity coefficient y+ as a function of cl’* for a 1:l electrolyte when a = 4.25 A, E = 78.7, and T = 298 K. The two upper curves are results from coupling parameter integration in the DH theory (eq 95) (- - -), labeled “DH 2”, and modified DH theory (eq 94) (- - -), labeled “MDH’. The data from GCMC simulations shown as (B) are taken from ref 23. The next two curves are from thermodynamic integration of the osmotic coefficient from the DH theory (eq 97) (...),labeled “DH(@)”,and modified DH theory (eq 96) (- - -), labeled “MDH(q5)”. The prediction of the usual DH expression for ln(y*),Le. the second term of eq 95, is shown as (- -), labeled “DH l”, while the DH limiting law is shown as (-), labeled “DH lim”.
-
are now obtained by inserting the pair distribution from eq 77 into eqs 84 and 86. The resulting integrals must be calculated numerically. For a symmetric electrolyte we have, for example,
where the second term on the right hand side is the contact term. The mean activity coefficient is then obtained from a numerical integration using eq 89. The excess intemal energy for 1:1 electrolytes in the modified DHX theory is virtually identical to that of the modified DH theory in Figure 7, so it is sufficient to use the linear theory (results from the nonlinear theory are not shown here). On the other hand, the nonlinear theory is essential to obtain good values of 4 for 1:l electrolytes, as mentioned above. In fact, the ordinary DHX theory gives very good values of 4 for 1:l electrolytes, as shown by Card and V a l l e a ~ . Our ~ modified DHX theory gives equally good results below 0.5 M. At higher concentrations it becomes unreliable and then breaks down since the potential of mean force only contains one exponentially decaying term (cf. discussion above for the linear case). The ordinary DHX approximation actually gives good values of 4 up to about 1 M concentration in the 1:l case: but this must be due to canceling errors since the DHX distribution function has the wrong functional form there. Since the values of 4 are very good in the ordinary DHX approximation for 1:1 electrolytes, the mean activity coefficient obtained from eq 89 is also very good, as shown by Card and V a l l e a ~ .The ~ same applies to the modified DHX theory below 0.5 M. Also in this case, the good results of the ordinary DHX theory up to about 1 M concentration4 are most likely due to canceling errors. We conclude that for 1 :1 electrolytes our modified DH and DHX theories give clear improvements in the values of the intemal energy compared to the ordinary DH and DHX theories. The osmotic and activity coefficients are predicted about equally well in the two sets of theories below 0.5 M. Since the ordinary DH and, in particular, the ordinary DHX theory give thermo-
Modified Debye-Huckel Approximation
J. Phys. Chem., Vol. 99,No. 25, 1995 10405 namely, to show how dressed ion theory and the concept of effective charge can be used to give a very simple, but efficient, improvement of the DH and DHX theories.
0.0 -0.5 -1.0
6. Summary
% l
-1.5 61
-2.0 -2.5 -3.0 0.0
0.4
0.2
0.6
0.8
cV2I MV2 Figure 9. Uex/NkBTas a function of cl’* for a 2:2 electrolyte of ionic diameter a = 4.2 8, when E = 78.6 and T = 298 K. Results from DHX theory (- * -), labeled “DHX’, from the modified DHX theory (-), labeled “MDHX’, and from GCMC simulations (B) are shown. The latter are taken from ref 23. I
I
-3
I
I
0.6
0.8
I 0.0
0.2
0.4 / Mu2
c112
Figure 10. Logarithm of the mean activity coefficient y* as a function of cl’* for a 2:2 electrolyte of ionic diameter a = 4.2 8, when E = 78.6 and T = 298 K. Results from DHX theory (- -), labeled “DHX’, from the modified DHX theory (-), labeled “MDHX”, and from GCMC simulations (B) are shown. The simulation data are taken from ref 23.
-
dynamic data for 1:1 electrolytes that are already quite good, this is to be expected. Nevertheless, since our modified theories are very simple, the improvements are worthwhile. For 2:2 electrolytes the linear theories and the ordinary DHX theory are clearly in~ufficient.~~ The results for the excess intemal energy from our modified DHX theory is presented in Figure 9. They are compared to GCMC simulation results taken from ref 23 and results from the ordinary DHX theory. Compared to the latter, we see that the modified DHX theory gives a clear improvement. The mean activity coefficient for 2:2 electrolytes is plotted in Figure 10, and we see that the modified DHX theory gives very good agreement with the simulation results up to about 0.25 M. Then it breaks down for the same reasons as in the 1:l case. Comparisons between Uex/NkBT and In p i results for 2:2 electrolytes from GCMC calculations and from several theories that are more complex than our modified DHX theory are given in ref 24. Inspection shows that compared to our theory, only the HNC and “ORPA B2” (a variant of the optimized random phase approximation) theories give equally good or better values of In &, while these two theories and the MPB theory give equally good or better values of Uex/NkBT (all this for concentrations less than 0.25 M). We can therefore conclude that our simple approach gives remarkably good results for 2:2 electrolytes, which are not easy to treat. There exist many other theories for electrolytes to compare with, but the comparisons we have made here suffice for the objectives of this paper,
+
We consider the distribution of ions around an arbitrary ion, the central ion, in a bulk electrolyte solution. To extend the Debye-Huckel theory, it is fruitful to regard the charge distribution in the ion cloud around the central ion as (essentially) a polarization charge. In the general case, the total charge that results from the polarization response in an electrolyte due to the electric field from some source charge equals
+
total charge = source charge polarization charge In our case, the central ion is the source charge. The average electrostatic potential is the potential that results from the total charge. This potential and the polarization charge distribution are self-consistent since the total polarization results from the field from both the source charge and the polarization charge. Dressed ion theory provides a simple, but formally exact generalization of the polarization response used in the DH theory (eqs 30b and 30a in DH correspond to eqs 34a and 34b in DIT, respectively). Since the DH theory is a linear approximation, it only contains a linear response. In the general case the polarization response is not linear, but we split it into a linear and a nonlinear part. The former part, &r), is expressed in DIT in terms of the average electrostatic potential by means of a self-consistent, linear response function of polarization, @(r) (eq 32). This response function (eq 25) is closely related to the dressed ion charge eo(r),which can be seen as the effective diffuse charge that interacts with the average potential (eq 32). If the nonlinear part of the polarization response together with the source charge is taken as an effective source charge, we can write
+
total charge = effective source charge polarization charge (linear response)
This has the same structure as in the DH theory, except that an effective source charge is used. The analogy to the DH theory is, however, very close since this effective charge is identical to the dressed ion charge (cf. eq 34). Thus, Qo(r) is both the effective source charge of the average potential and the effective charge that interacts with the average potential, two perfectly symmetric roles. We conclude that Q0(r) in the exact dressed ion theory takes the roles that the bare ionic charge has in the DH approximation. The interaction energy between the dressed ion charge and the average potential (we denote this energy -h#(r)/P; see eq 36) is the leading term in the potential of average force w(r) at low concentrations. By neglecting all other contributions to w(r) (eq 67), we obtain a generalization of the PoissonBoltzmann theory (eq 54) that-in contrast to the ordinary PB theory-has symmetric pair distribution functions g, = gJ,(eqs 53 and 68). Here the dressed ion charge replaces the bare ionic charge used in the ordinary PB theory. In contrast to DIT, this symmetric PB theory is an approximation. If we, in addition to this approximation, assume that h#(r) has the same behavior at all r values as it has at large r (eqs 69, 75, and 76), we obtain the modification of the nonlinear DH theory (modified DHX theory) proposed in this work. In the linear theory we furthermore linearize the Boltzmann factor in g(r) (the exponential function in eq 68) and obtain eq 70. The resulting modified DH and DHX theories differ from the usual DH and DHX theories by using effective point charges instead of bare ionic charges and by enforcing the exact relationship
10406 J. Phys. Chem., Vol. 99, No. 25, 1995 (41) for the decay parameter K . This leads to K/KD and q*/q values as functions of concentration for 1:1 and 2:2 electrolytes at low to moderate concentrations in quite good agreement with those from the HNC approximation (note that q*/q = ( K / K D ) ~ for symmetric electrolytes). For 2:2 electrolytes it is, of course, necessary to use the nonlinear theory, and then the characteristic dip of K/KD and q*/q below the value 1 at low concentrationsI2 is fairly well reproduced. The ordinary DH and DHX theories as well as the ordinary PB theory have, of course, K = KD throughout. In our modified DH theory K/KD as a function of K D is~ given as a solution of a simple transcendental equation (eq 65b). Since it is independent of the second dimensionless parameter that defines the system, / 3 q 2 / ( 4 n e ~ (which ~) is denoted p*), it has the same erroneous behavior at low concentrations or small ion sizes as, for example, in the MSA and LMPB theories. This is not the case in our modified DHX theory, where K/KD is a function of both KDa and p*. The main mechanisms for the behavior of q* (and hence of K ) that are includedfor all ions in our modified DHX treatment are (i) the exclusion of ion cloud charge from the hard core region, which leads to the buildup of an appreciable charge density outside the core, and (ii) the nonlinear polarization response from the action of the average potential. The latter effect is particularly important when p* is large or K D is~ small, and it gives a negative contribution to q*/q. The effect (i) gives a positive contribution to q*/q and is a dominant effect for 1:l electrolytes except at very low concentrations, where (ii) dominates. In the linear theory only (i) is included. In our modified DH theory the dressed ion charge Bo(r) is entirely located inside the core region of each ion (eq 74), while in our modified DHX theory there is a short-range contribution outside the core as well (eq 81). The role of the charge contribution inside the core is to cancel the charge that would be placed there by the average electrostatic potential if the hard core exclusion were absent. (This is an indirect way to enforce the hard core exclusion that gives the positive contribution to q*/q mentioned above.) The contribution to eo(r)outside the core in the modified DHX theory is the nonlinear polarization response due to the average potential. Note that in the exact DIT there are other contributions to so(r) which affect the behavior of K/KD and q*/q (see refs 11 and 12), but the effects (i) and (ii) are major contributions at low to moderate concentrations. The excess internal energy, osmotic coefficient, and activity coefficient calculated in the modified DH theory-or, when appropriate, the modified DHX theory-agree well with those from Monte Carlo simulations at low to moderate concentrations. In the modified DH theory there are explicit formulas for the thermodynamic quantities in terms of K , closely analogous to the formulas in the ordinary DH theory. Acknowledgment. This work received financial support
from the Swedish Natural Science Research Council. Appendix A
Kjellander
+
of species i at r2 changes from ni(rz) to, say, ni(r2) dni(r2). If Bvj(r1) for a l l j is sufficiently small everywhere, the change in density dni(r2) is linearly related to dvj(r1) (nonlinear contributions can be neglected when Bvj(r1) 0 for all rl). The Yvon equation2 states that
-
dn,(r,) = -pn,(r2) dvi(r2) -
For a bulk fluid we have vj(r1) = 0 everywhere for all j , the density is independent of the coordinate, nl(r) = nl = constant, and the correlation functions are spherically symmetric, hlm(r,r’)= hlm(lr- r’l). Consider a bulk electrolyte where we tum on an extemal electrostatic potential dYext(r)which is small everywhere. Then we have 6vj(r) = a BYext@), and by applying eq A1 we obtain eq 15. We next turn to the relationship between h,(r) and hUo(r). The direct correlation function cl,(r)is defined by the OmsteinZemike (OZ) equation,*
from which q ( r ) can be determined once h&) is known. For large r values q ( r ) decays to 0 proportionally to the pair potential uij(r); we have cu(r) -puij(r). For ions, only the electrostatic potential contributes at large r, and we have
-
Therefore, the function c,“(r) defined by
is a more short-ranged function than cv(r). The function huo(r) is defined” from cdo(r)via the OZ equation,
and can accordingly be determined from hij(r) via eqs A5, A4, and A2. Finally, we shall look more closely at the potential of mean force, wv(r). The bridge function ev(r) can be defined from
+
+
w&r) = -/3uij(r) h&r) - cij(r) e&r) (A6) and it is possible to write it as an infinite sum of complicated, multidimensional integrals (bridge diagrams) that only contain products of the total correlation functions and the densities.25 If we insert eqs A4 and 38, we can write eq A6 as wij(r)= - p u T ( r )
+ h,#(r) + h:(r)
- ci:(r)
+ e&r) (A7)
where uycore is the hard core potential which is 0 for r Comparing with eq 51, we can identify
+
f j ( r ) = hijo(r) - cij0(r) ev(r),for r L a
In this appendix we shall summarize some facts referred to in the main text that are somewhat more technical in nature than the rest of this paper. We start with a result from linear response theory. Consider a fluid where the particles interact with an external object. Let v&) denote the interaction potential between the object and each particle of species I. The density distribution of species 1 in the fluid is nr(r), and the total correlation function is hm(r,r’). If the extemal potential for speciesj at rl is changed from vj(r~)to v,(r~) dvj(r1) (for all j and all rl). the density
+
-
2
a.
(~8)
-
When the ion concentration goes to 0, it follows from eq A5 that hio(r) - cio(r) 0 for all r, and we also have ev(r) 0, so we can conclude thatfij(r) 0 for all r. In the HNC approximation we set ev(r) = 0, so
-
f , ( r ) = hi;(r) - c:(r)(HNC)
(-49) and we see that the approximation in eqs 53 and 54, where fil(r) = 0 for all r, is somewhat simpler than the HNC approximation.
Modified Debye-Hiickel Approximation
Appendix B In this appendix we shall prove eq 63 in the modified DH theory. Using eq 60, we can write eq 57 as K2q.*e-Kr I @,(r) = - ,rTa 4nr We assume that the total correlation function has the mathematical form h&r) = a , a F ( r ) (B2) with some constants ai and aj and some function F(r). It follows from eqs 1 and B2 that
This should be identified with eq B 1. We may therefore-without loss of generality-take F(r) = -AB exp(-~r)/[4m~or],where A is some constant, and identify
By multiplying both sides of this equation by @ n,, summing over j , and inserting eq 41, we obtain after simplification
Equations B4 and B5 give and eq B2 now yields eq 63.
References and Notes (1) Debye, P.; Hiickel, E. Phys. Z. 1923, 24, 185.
J. Phys. Chem., Vol. 99, No. 25, 1995 10407 (2) Hansen, J. P.; McDonald, I. R. Theory of Simple Liquids, 2nd ed.; Academic Press: London, 1986. (3) Pitzer, K. S. Arc. Chem. Res. 1977, 10, 371. (4) Card, D. N.; Valleau, J. P. J. Chem. Phys. 1970, 52, 6232. ( 5 ) See, for example, ref 3 and the following: Conway, B. E. In Comprehensive Treatise of Electrochemistry; Thermodynamic and Transport Properties of Aqueous and Molten Electrolytes; Conway, B. E., Bockris, J. O'M., Yeager, E., Eds.; Plenum Press: New York, 1983; Vol. 5 , p 111, and references therein. (6) Outhwaite, C. W. In Statistical Mechanics: A Specialist Periodical Report; The Chemical Society: London, 1975; Vol. 2, p 188. (7) Outhwaite, C. W.; Molero, M.; Bhuiyan, L. B. J. Chem. Soc., Faraday. Trans. 1993, 89, 1315. (8) Nordholm, S.Aust. J. Chem. 1984, 37, 1. (9) Nordholm, S. Chem. Phys. Lett. 1984, 105, 302. Penfold, R.; Nordholm, S.;Jonsson, B.; Woodward, C. E. J. Chem. Phys. 1991, 95, 2048. (10) Kjellander, R.; Mitchell, D. J. Chem. Phys. Lett. 1992, 200, 76. (11) Kjellander, R.; Mitchell, D. J. J. Chem. Phys. 1994, 101, 603. (12) Ennis, J.; Kjellander, R.; Mitchell, D. J. J. Chem. Phys. 1995, 102, 975. (13) Attard, P. Phys. Rev. E 1993, 48, 3604. (14) Leote de Carvalho, R. J. F.; Evans, R. Mol. Phys. 1994, 83, 619. (15) In refs 10 and 11 is denoted -a,and in ref 10 we use the notation hS for ho and ehntfor eo. (16) Kirkwood, J. G. Chem. Rev. 1936.19.275, Kirkwood, J. G.; Poirier, J. C. J. Phys. Chem. 1954,58, 591. (17) FBat, G.; Levine, S. J. Chem. Sac., Faraday. Trans. 2 1977, 73, 1345. (18) Schmidt, A. B. Soviet Electrochem. 1984, 20, 767. (19) Outhwaite, C. W. J. Chem. SOC., Faradny. Trans. 1991, 87, 3227. (20) Outhwaite, C. W. Chem. Phys. Lett. 1970, 5, 77. (21) Waismann, E.; Lebowitz, J. L. J. Chem. Phys. 1!?72,56,3086,3093. (22) McQuarrie, D. A. Statistical Mechanics; Harper and Row: New York, 1976. (23) Valleau, J. P.; Cohen, L. K. J. Chem. Phys. 1980, 72, 5935. (24) Valleau, J. P.; Cohen, L. K.; Card, D. N. J. Chem. Phys. 1980, 72, 5942. (25) Stell, G. In The Equilibrium Theory of Classical Fluids; Frisch, H. L., Lebowitz, J. L., Eds.; W. A. Benjamin: New York, 1964; pp 11-171. JF950023D
