Volume
127, number
INTERACTION Roland
CHEMICAL
4
OF CHARGED
KJELLANDER
PHYSICS
SURFACES
and Stjepan
LETTERS
IN ELECTROLYTE
20June1986
SOLUTIONS
MARCELJA
Department of Applied Marhematlcs, Research School of Physical Sciences, Austrahan National University, Canberra, ACT 2601, Australia Received
11 February
1986: m final form 2 April 1986
We report a restricted primitive model calculation of the double-layer interaction between two uniformly charged surfaces immersed in an electrolyte solution, where the anisotropic hypernetted chain approximation is utilized for the pair correlations. The strong, attractive pressure contribution resulting from ion-ion correlations is very similar to that found earlier for double layers with only counterions present. Consequently, the Poisson-Boltzmann equation, which neglects ion-ion correlations, significantly overestimates the double-layer repulsion in most situations. Some ion concentration profiles are presented, and for large separations between the surfaces the profiles close to each surface agree with grand canonical simulation results for a single wall.
1. Introduction The commonly accepted description of electric double layers is still the classical picture based on the non-linear Poisson-Boltzmann (PB) equation. The equation is derived for a model where the solvent is a dielectric continuum and the ions are point charges, and it entirely neglects the ion-ion correlations. However, it has recently been shown via Monte Carlo simulations [l] and the anisotropic hypernetted chain (HNC) approximation [2] that, for the same model, the solutions of the PB equation are quantitatively, and sometimes even qualitatively, incorrect in describing the interaction between charged surfaces. The ionion correlations give rise to a large, attractive pressure contribution of electrostatic origin, which, for some conditions, may even dominate in the total pressure. The work in refs. [ 1,2] has been restricted to model systems where point-charge counterions are the only ionic species present in the solvent. Since the ions repel each other electrostatically, the point-charge description is quite accurate provided the real ionic radius is not too large. For fairly large ions, a repulsive pressure contribution due to the core-core interactions becomes important [3]. This contribution compensates to some extent the attractive, electrostatic pressure term. For large ion radii it eventually dom402
inates over the attractive term. In the presence of an electrolyte, it is essential to have a finite ionic size. The ionic core repulsion prevents the positive and negative ions from collapsing on each other. We have chosen to describe the core potential as a hard-sphere interaction, and our system corresponds to a restricted primitive model electrolyte between uniformly charged surfaces. (For bulk solutions, hard-core interactions and short-range, soft-core interactions have been shown to lead to very similar results [4] .) The results reported below are, to our knowledge, the first calculations of the interaction between charged surfaces in a primitive model electrolyte solution, where the level of approximation is sufficiently high to describe accurately the inhomogeneous phase between the surfaces. We find that most of the important features observed in the counterion-only systems are retained in the presence of the electrolyte, and that the interaction calculated via the PB equation in general does not agree with the accurate solutions for the same model system.
2. Method Our calculation
is based on a strict equivalence
of
0 009-26 14/86/$ 0350 0 Elsevier Science Publishers B.V. (Norht-Holland Physics Publishing Division)
Volume 127. number 4
an inhomogeneous fluid between two planar walls with a multispecies two-dimensional homogeneous fluid [S] . In case of ionic fluids, we use the anisotropic HNC approximation to close the integral equations for the paircorrelation functions. The HNC approximation is very accurate in treating pair correlations due to Coulomb interactions. It is less accurate for the hard-core interactions, but results for bulk electrolyte solutions [6] show that the accuracy is nevertheless still very good at densities typical of electric double layers. Note that the anisotropic HNC approximation is formulated on the two-particle distribution level. It is a more sophisticated theory than the various double-layer theories that utilize bulk ion-ion correlation functions (usually the direct correlation function) to calculate the one-particle distributions in the inhomogeneous phase. In our approach we explicitly calculate both the anisotropic ion-ion correlation functions (in presence of the walls) and the concentration profiles [S]. In order to equilibrate the inhomogeneous ionic fluid between the surfaces with a bulk electrolyte and to evaluate the net pressure between the surfaces, we need bulk values of the chemical potential and the osmotic coefficient. The values for typical primitive model electrolytes in the NIX approximation are available in the literature [7]. Once the bulk chemical potential is set, our method very easily adjusts the ion concentrations between the surfaces until the chemical potential at each position is equal to that of the bulk [5]. The ionic number density between the surfaces equals P”(Z)
=
exp@,) q?,(z)
20June1986
CHEMICAL PHYSICS LETTERS
(1) ’
where the coordinate z is counted perpendicularly to the surfaces and where V, is the chemical potential for species v (index v here assumes the values + or -), /3 = l/k,T, k, is Boltzmann’s constant, T the temperature, LJ, the thermal wavelength and y,(z) the activity coefficient at location z. A very great advantage with the HNC approximation is that the local activity coefficient can be explicitly evaluated from the correlation functions without performing any numerical coupling constant integration [5,8]. Therefore, eq. (1) can easily be used to obtain a self-consistent density profile once 12, is known [5]. In the present case, the bulk chemi-
cal potential
cc+ bulk is given, and we accordingly know (a symmetric electrolyte is as:oL++II_)=P,bulk sumed). This gives P+ = fi!” t $ and cc_ = pplk J/ with an undetermined constant $, and we obtain
P,(Z)
expC4ErPlk) =_p exp@/-Pk)(2) =ew(@lCI) A3Yy,(Z) .A3Y,(4 ’
where X = exp@Jl) is as yet unknown and A = A+ = A _. The electroneutrality condition for the system between the surfaces gives the following simple, quadratic ic equation for X:
( J’ 4+ X
-6 y+(z)
+_!_ x
’ 4_
h
J Y_(Z) ) -6
exp@~tT A3
t20=0, (3)
where q+ and q_ are the ionic charges, u is the surface charge and 26 is the surface separation. (For an asymmetric electrolyte, we would obtain a similar, higherorder equation.) Eqs. (2) and (3) together with the HNC equations for the anisotropic pair correlation functions form a complete set of equations that can be solved self-consistently [5]. A technical point, worth mentioning, is that the long-range tails of the various kinds of pair-correlation functions in the real space and in Fourier space, as induced by the Coulomb and the hard-core interactions, are subtracted and then treated analytically in our computer implementation of the scheme. This is important in order to achieve a sufficient numerical accuracy. The results in this paper were obtained on a VAX 1 l/750. Up to 122 discrete layers between the surfaces and 100 to 150 mesh points in the lateral direction were used.
3. Results and discussion We have tested the performance of the anisotropic HNC scheme and the accuracy of our computer program by comparing our results for two surfaces at a fairly large separation with the grand canonical Monte Carlo data for a single surface as calculated by Torrie and Valleau [9]. Fig. 1A shows the ion density profiles when the bulk phase is a 1.O M 1: 1 electrolyte 403
Volume 127. number 4
CHEMICAL PHYSICS LETTERS
20June1986
A
‘a DISTANCE 6,
Fig. 1. (A) Comparison between ion density profiles as obtained by the grand canonical Monte Carlo technique (e counterions, A coions; taken from ref. (91) and by the anisotropic HNC approximation (full curves). The MC results are for a single charged wall (U = 0.620 C/m*) immersed in a 1 .O M 1: 1 electrolyte solution, while the HNC results are for two such walls at fairly large separation (40 A), where the intervening solution is in equilibrium with a bulk solution of the same composition as in the MC calculations. The ionic diameters are 4.25 A. Only the left halves of the profiles are shown in the two-wall case; the coordinate 0 then shows the midpoint between the walls. For one wall, the electrolyte solution continues to the right and turns into a bulk phase. (The distance scale is chosen to suit two walls.) The point of closest approach by the ions to the left surface is located at the coordinate -20 A in both cases. (B) The results given by the anisotropic HNC approximation for the same two-wall system as in (A), but for a separation of 15 A. The upper curve shows the counterion profile while the lower shows the coion profile. For all cases in this figure the temperature is 298 K and the dielectric constant of the solvent is 78.5. The bulk chemical potential is taken from the MC calculations.
solution and the surface charge is 0.620 C/m2. Although such high surface charges are seldom encountered in practice, we have chosen that example because its theoretical description is difficult. Our method has successfully reproduced the second density peak at about 4 A from the surface. The peak is a consequence of the crowding of the ions close to the surface and it indicates an ion layering (cf. refs. [9,10]). Fig. 1B shows the profile for the two-surface case at a smaller Fig. 2. The counterion (A) and the coion (B) density profiles b for the system between two charged surfaces (D = 0.224 C/m*) in equilibrium with a 0.5 M 1 :l electrolyte solution. The full curves show the HNC results, while the dotted curves show the predictions from the non-linear Poisson-Boltzmann theory. The surface separation is 15 A, and the ionic radii are 2.3 A. The temperature is chosen 298.15 K and the dielectric constant of the solvent is 78.358. Note that the ordinate scales of the two figures are very different.
404
DISTANCE
tk
20 June 1986
CHEMICAL PHYSICS LETTERS
Volume 127, number 4
separation, where the two double layers overlap substantially. The structure of the profile is then increased. The density profiles at a lower surface charge density and lower electrolyte concentration are shown in fig. 2. The counterion profile is very similar to the corresponding PB profile, but one can seen that the ionion correlations make the counterions come closer to each surface. This is similar to the situation for the case with only counterions present [ Ill. On the other hand, the coion profile deviates considerably from the corresponding PB profile, as shown in fig. 2B. These ions tend strongly to avoid the surfaces. Perhaps the most interesting of our new results are those for the interaction between two charged surfaces immersed in an electrolyte solution. The surface-surface interaction (the perpendicular component of the pressure tensor) is calculated as a sum of distinct physical contributions:
:
1
.. :.. *. :
..*.
‘5. --.. X. L. . .. . **.* . . ..
where Pkin = kBT C P,(O) 3 ”
pel=mvq,j dzPv(z)j!di’P,.(r) -6
&F’u’(r,z, z’)
xsd’ vu az P core = k,T
c lJ,v’
2n j
da
h,,4, z, z’) p,(Z)
j
9
x(Z-z’)g”,~([D2-(z-z’)2]1’2,z,z’), Pb,,lk = kBTqbbulkc ”
-..
a...
-..
-..
....
(8,
Fig. 3. The net pressure between two charged walls (u = 0.224 C/mZ) immersed in a 0.5 M 1 :l electrolyte solution as a function of surface separation (lower curves). As a comparison, the pressure between the walls when only counterions are present is also shown (upper curves). The full curves show the HNC results, and the dotted curves show the PB predictions (log scale). The surface separation is defined as the distance between the points of closest approach of the ions to the surfaces. (The remaining parameter values are given in the caption of fig. 2.) For surface separations below about 2.7 A the HNC pressure (cf. fig. 4) becomes smaller than the PB result.
dZ’ /+(z’)
-s
0
.. . .
:
SEPARARATION
0
-.*. -1..
\\
’
*-..
(4)
+pel + ‘core - Pbulk 3
‘tot = pti
‘..*
t* :
pE”lk .
u”?’ is the normal Coulomb interaction, h,,, =gVV, -Y ,gvv* 1s . the pair distribution function, r = (x,y), D is the hard-core ion diameter and Gbulk is the bulk osmotic coefficient (taken from ref. [7]). The integrals in pore are only to be taken for the values of z and z’ where the argument of the square root is positive. An example for a 0.5 M 1: 1 electrolyte and a surface charge of 0.224 C/m2 is shown in fig. 3, where it is compared with results for the similar system with only counterions between the surfaces. The corresponding PB pressures are also shown. With our current
numerical accuracy, in this example we could only compute the HNC pressure for relatively small surface separations in presence of electrolyte. In the asymptotic regime at larger separations, the different terms in eq. (4) (each of which is relatively large) will almost completely cancel each other. However, when the electrolyte concentration is decreased, the pressure can be accurately calculated for larger surface separations because it does not decay as rapidly. The computed results show that at small surfacesurface separations the interactions obtained when counterions are the only species in the solution do not become substantially changed by the addition of the electrolyte. This is not unexpected, because under such conditions the total number of coions in the double layer is small. The decrease in the net force when 405
Volume 127, number 4
CHEMICAL PHYSICS LETTERS
20 June 1986
I 3-
0
1
6
2
7
8
0
5
SE&AR*4TIoN:h
15 SERm&N
tk
Fig. 4. The different physical components (defined in eq. (4)) of the net pressure as calculated by the HNC method for the same case as shown in figs. 2 and 3. The regions of small or larger surface separations are shown in (A) and (B) respectively. The bulk pressure to be subtracted is Pb&R T = 1.026 M.
the electrolyte concentration is raised, is mainly due to the increase in Pbulk . It is particularly interesting that the addition of the electrolyte has a relatively small effect on the electrostatic attraction term. For example, at the surface separation of 15 A, Pel increases (in absolute value) by about 5% upon addition of electrolyte, while at the same time the total repulsive pressure is decreased by a factor of four. The contributions to P,l from the coion-counterion interactions are to some extent compensated by the coion-coion and counterioncounterion interactions. On the other hand, the added substantially, mainly through electrolyte changesp,,, the coion-counterion interactions (no compensation is possible because all contributions are positive). At the surface separation 1.5 A, the increase inP,,,, is 70%. The corresponding numbers at the separation 20 A are 40% for Peland more than a factor of 4 for Pare. The different contributions to the pressure between the surfaces for the example discussed above are shown in fig. 4. It is remarkable that in the regime where electrostatic and core contributions are as large as the kinetic term, the PB result, which only considers a ki406
netic pressure, is reasonably good. However, in the case of a divalent electrolyte, our preliminary results indicate that the electrostatic attraction becomes dominant, and just as in the case where only counterions are present [l-3] the double-layer repulsion turns into an attraction. Acknowledgement As in our previous numerical work we have been kindly supported by David Sandilands and David Smith from the Computer Unit at the Research School of Biological Sciences, ANU. We also thank Phil Attard for programming the PB equation routines.
References Ill L. Guldbrand, B. Jiinsson, H. Wennerstrom and P. Lime, J. Chem. Phys. 80 (1984) 2221.
121 R. Kjellander and S. Mar?elja, Chem. Phys. Letters 112 (1984)49;
114 (1985) 124E.
131 R. Kjellander and S. Marklja, J. Phys. Chem., to be published.
Volume 127, number 4
CHEMICAL PHYSICS LETTERS
[4] P.J. Rossky, J.B. Dudowicz, B.L. Tembe and H.L. Friedman, J. Chem. Phys. 73 (1980) 3372. [S] R. Kjellander and S. Mari?elja, J. Chem. Phys. 82 (1985) 2122. [6] J.C. Rasaiah, D.N. Card and J.P. Valleau, J. Chem. Phys. 56 (1972) 248; J.P. Valleau, L.K. Cohen and D.N. Card, J. Chem. Phys. 72 (1980) 5942. [7] J.C. Rasaiah and H.L. Friedman, J. Chem. Phys. 48 (1968) 2742. [8] J.P. Hansen, G.M. Torrie and P. Viellefosse, Phys. Rev. Al6 (1977) 2153.
20 June 1986
[9] G.M. Torrie and J.P. Valleau, J. Chem. Phys. 73 (1980) 5807; S.L. Carnie and G.M. Torrie, Advan. Chem. Phys. 56 (1984) 141. [lo] P. Nielaba and F. Forstman, Chem. Phys. Letters 117 (1985) 46; C. Caccamo, G. Pizzimenti and L. Blum, Phys. Chem. Liquids 14 (1985) 311. [ 111 B. Jtinsson, H. Wennerstrom and B. Halle, J. Phys. Chem. 84 (1980) 2179.
407