4210
Langmuir 2002, 18, 4210-4219
Comparison of the Thermal and Concentration Criteria for Determining the Effective Charge of a Spherical Colloidal Particle in a Salt-Free Suspension Arup K. Mukherjee,† K. S. Schmitz,‡ and L. B. Bhuiyan*,† Laboratory of Theoretical Physics, Department of Physics, University of Puerto Rico, Box 23343, UPR Station, San Juan, Puerto Rico 00931-3343, and Department of Chemistry, University of Missouri-Kansas City, Kansas City, Missouri 64110 Received November 1, 2001. In Final Form: February 21, 2002 The techniques of canonical Monte Carlo simulations are employed in a comparative study of two recently proposed criteria for the determination of the renormalized or effective charge of a spherical macroion in a charge-stabilized colloidal suspension. The effective charge is determined using (i) the thermal criterion, where the effective charge consists of the intrinsic macroion charge reduced by the integrated thermodynamically bound counterion charge, the bound counterions being defined as those counterions that are within and up to the radial distance (from the macroion center) at which the reduced mean electrostatic potential energy equals the thermal energy, and (ii) the concentration criterion, where the bound counterions are those that lie within the distance at which the local concentration of the counterions equals its average value. It is found that the former criterion applies to the entire range of bare macroion charge/macroion radius ratios investigated, whereas the latter criterion is more suited to large values of this ratio. The observation that the effective charge/bare charge ratio (Zeff/Zm) depends only on the bare charge/macroion radius ratio (Zm/rm) suggests that all macroions with the same numerical value of the Zm/rm ratio are “equivalent systems” and that one can therefore study highly charged colloids by using a more manageable set of parameters in a simulation. The calculations are extended to divalent counterion systems, and comparisons are also made with other criteria for determining the effective macroion charge in the literature.
1. Introduction More than half a century after it was introduced, the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory for the pair potential still enjoys today the classic status of being the most widely used theory for the interpretation of the stability of colloids.1 The DLVO potential is composed of two parts, an attractive van der Waals term (UvdW) that is significant only at short distances between two particles and a long-range screened Coulomb repulsive part (UC). Although the form of the potential is available for planar and spherical geometries, the focus of the present study is on spherical particles. Since the inception of the DLVO potential, there have been many variations on the form of the pair potential. The general form of the pair potential for two particles i and j of common radius rm and charges (with sign) Zi|e| and Zj|e| (|e| being the magnitude of the electronic charge) can be written as
βUij ) βUC + βUvdW
[
rm2 rm2 exp(-κr) βAH ) λBZiZjBiBj + + r 3 r2 - 4r 2 r2 m
(
)]
r2 - 4rm2 1 loge 2 r2
(1)
where β ) 1/kBT; kB is the Boltzmann constant; T is the absolute temperature; λB ) βe2/(4π0r) is the Bjerrum length; 0 is the vacuum permittivity; r is the relative * To whom correspondence should be addressed. † University of Puerto Rico. ‡ University of Missouri-Kansas City. (1) Verwey, E. J. W.; Overbeek, J. Th. G. Theory of the Stability of Lyophobic Colloids; Elsevier: Amsterdam, 1948.
permittivity (bulk dielectric constant) of the solvent; r is the center-to-center distance between the macroions; AH is Hamaker’s constant, which lies in the range 10-12 J < AH < 10-19 J; and the Bm are model-dependent parameters for the mth macroion. The screening parameter κ is given in generalized form by
κ2 ) 4πλB(
∑j njZj2 + Zc2nc) ) κadd2 + κc2
(2)
where Zj and nj are the valency (with sign) and number density, respectively, of the jth species of added electrolyte and nc and Zc are the number density and valency (with sign), respectively, of the counterions that originated from the macroions. We have distinguished contribution of the the added electrolyte to the screening parameter (κadd) from that of the counterions (κc) because the former can be controlled independent of the macroions and the latter represents the lower limit of κ, as it depends solely on the macroion concentration. The DLVO potential is recovered by the substitutions Bm ) exp(κrm)/(1 + κrm) and κ ) κadd. That is, the DLVO potential is derived as a limiting law expression for the interaction between two isolated macroions in the presence of “excess” added electrolyte. To our knowledge, the first paper to point out that the pair potential depends on other macroions present in the system is by Beresford-Smith and co-workers (BSCM model).2 Using the Jellum approximation (uniform colloid distribution), they obtained a volume fraction (φm) dependence for the factor Bm, namely, Bm ) (1 + φm) exp(κrm)/(1 + κrm), for the linearized PB solution. Belloni3 obtained a solution to the BSCM model without the Jellum (2) Beresford-Smith, B.; Chan, D. Y. C.; Mitchell, D. J. J. Colloid Interface Sci. 1985, 105, 216. (3) Belloni, L. J. Chem. Phys. 1986, 85, 519.
10.1021/la011625t CCC: $22.00 © 2002 American Chemical Society Published on Web 04/24/2002
Effective Charge of a Particle in a Salt-Free Suspension
approximation, with the result that Bm ) cosh(κrm) + [κrm cosh(κrm) - sinh(κrm)][(z/κrm3) - (γ/κrm)], where z ) 3φm/(1 - φm) and γ is a function of κ, Zm, and rm, whose exact form is not important in the present discussion. The Yukawa potential is obtained for Bm ) 1. Our interest in the present study lies only in the Coulomb part of the pair potential; hence, we omit any further discussions of the van der Waals attraction term. It is evident from the early work on the pair potential that theories based on the linearized Poisson-Boltzmann (PB) solutions give rise to a screened Coulomb contribution of the form of UC in eq 1. We note further that the validity of the first term in eq 1 was examined by Stigter,4 who showed that the screened Coulomb form is valid only in the outer part of the overlapping double layers. More recently, Stevens, Falk, and Robbins5 emphasized that the Yukawa form of the pair potential and nonlinear solutions to the PB equation might be valid for the interstitial regions but break down in the region near the macroion surface. Nonetheless, the form of UC has routinely been applied to colloidal and polyelectrolyte systems in which there is considerable overlap of the ion clouds of adjacent particles. For example, Robbins, Kremer, and Grest6 reported molecular dynamics simulations of spheres interacting through the Yukawa potential in which the criterion set forth by Stigter4 was violated, i.e., βUij > 1. Similarly, Monte Carlo (MC) simulations such as those of Lobaskin and Linse;7 Linse and Lobaskin;8 and Lobaskin, Lyubartsev, and Linse9 violated the Stigter criterion. A general result in applying eq 1 to experimental data was that the fitting value of Zm is considerably smaller than the value obtained by chemical methods. Roberts and co-workers10 recently tabulated experimental data on the ratio Zeff/Zm for a wide range of values of Zm and rm. They showed that the ratio Zeff/Zm decreased monotonically as the ratio Zm/rm increased in value. Thus, the concept of an effective charge, Zeff, or charge renormalization was born with the following apparent objective, as noted by Ulander, Greberg, and Kjellander:11 “The key idea is to maintain the feature from above that the propagation of the electrostatic potential in the electrolyte is governed by the screened Coulomb potential φo(r).” That is, the introduction of a renormalized charge provides a justification for the use of eq 1 by the reduction of the magnitude of Zm. One approach to obtaining the effective charge is to determine an effective pair potential between macroions in a simulation and then attempt to fit this curve with the assumed form of a screened Coulomb interaction. Lobaskin, Lyubartsev, and Linse9 recently reported MC simulations on a highly asymmetric charge system where the macroion charge was -60 in the presence of monovalent or divalent counterions. They found that the effective charge in the presence of divalent ions was less than that for monovalent ions when the macroion-macroion pair distribution function was interpreted in terms of a DLVOtype potential. This latter quantity was obtained through (4) Stigter, D. Rec. Trav. Chim. 1954, 73, 593. (5) Stevens, M. J.; Falk, M. L.; Robbins, M. O. J. Chem. Phys. 1996, 104, 5209. (6) Robbins, M. O.; Kremer, K.; Grest, G. S. J. Chem. Phys. 1988, 88, 3286. (7) Lobaskin, V.; Linse, P. J. Chem. Phys. 1999, 111, 4300. (8) Linse, P.; Lobaskin, V. J. Chem. Phys. 2000, 112, 3917. (9) Lobaskin, V.; Lyubartsev, A.; Linse, P. Phys. Rev. E 2001, 63, 020401-1. (10) Roberts, J. M.; O’Dea, J. J.; Osteryoung, J. G. Anal. Chem. 1998, 70, 3667. (11) Ulander, J.; Greberg, H.; Kjellander, R. J. Chem. Phys. 2001, 115, 7144.
Langmuir, Vol. 18, No. 11, 2002 4211
an inverse MC procedure on a system containing only macroions. However, it is not clear how much emphasis was placed on the role of the counterions vis-a`-vis the effective charge in these MC studies. For instance, the redistribution of these counterions upon approach of two (or more) macroions is not explicitly considered in the calculation of the effective charge. A second approach is to estimate an effective charge by means of a charge renormalization from the distribution of counterions about a central macroion5,11-17 Alexander et al.12 proposed a theory based on the equilibrium set up for electrostatic interactions of a counterion with a renormalized surface charge and its entropy of mixing with the solution, viz., ZeffλB/rm ≈ -ln(φc), where φc is the volume fraction of the counterions and κ ) κc. They fitted the pressure obtained by adjusting the macroion charge in the Debye-Hu¨ckel theory (viz., the first term in eq 1) to the pressure given by the nonlinear PB theory. The procedure was repeated recently by Stevens et al.5 and Lobaskin and Linse,7 but they used Monte Carlo (MC) simulations instead of the PB theory. Gisler et al.13 obtained a value for the renormalized charge by matching the nonlinear and linear solutions to the PB equation at the computation cell boundary and then extrapolating back to the macroion surface. An important result of both Alexander et al.12 and Gisler et al.13 is that Zeff reaches a saturation value for large Zm. Belloni14 has criticized this matching method of Gisler et al. because the screened Coulomb form mimics only the long-distance behavior of the electrostatic coupling. Somewhat related is the work of Ulander, Greberg, and Kjellander,11 who were concerned primarily with planar geometry but included some discussion of spherical macroions. As indicated in the quote above, their objective was to preserve the screened Coulomb form of the electrical interactions. They treated both the macroions and the electrolyte ions as dressed ions with an effective charge. Rather than use the solvent dielectric constant, they employed the solution dielectric constant in their analysis of the propagation of the electric field. What resulted were two effective charges that decayed in accordance with two screening parameters. Implicit in the renormalization schemes is the notion that neutralizing counterions are found in the vicinity of the surface of the macroion. In this way the macroions are allowed to sample the entire free volume defined by the quantity V(1 - φm), where the volume fraction φm is calculated from the hard-sphere radius of the macroion, rm, for example, φm ∝ rm3. It is curious to note that, in their seminal work, Verwey and Overbeek1 anticipated charge renormalization in their citing of the work of Mu¨ller15 on his numerical solution to the Poisson-Boltzmann equation. Noting that the potential decays rapidly from the surface of the sphere, Verwey and Overbeek stated on page 401 “Hence this part of the potential curve will have rather a small influence upon the form of the curve in the entire diffuse layer, as it comprises only a small part of the total space curve”. The above charge renormalization schemes are directed at a reduction in the surface charge density. That is, the distribution function about the macroion is highly “collapsed” onto the surface of the macroion. This idea appears (12) Alexander, S.; Chaiken, P. M.; Grant, P.; Morales, G. J.; Pincus, P.; Hone, D. J. Chem. Phys. 1984, 80, 5776. (13) Gisler, T.; Schulz, S. F.; Borkovec, M.; Sticher, H.; Schurtenberger, P.; D’Aguanno, B.; Klein, R. J. Chem. Phys. 1994, 101, 9924. (14) L. Belloni, L. Colloids Surf. A 1998, 140, 227. (15) Mu¨ller, H. Z. Phys. 1927, 24, 185. (16) Linse, P. J. Chem. Phys. 2000, 113, 4359 (17) Shklovskii, B. I. Phys. Rev. E 1999, 60, 5802.
4212
Langmuir, Vol. 18, No. 11, 2002
Mukherjee et al.
to have been echoed in the recent MC simulations of a highly asymmetric, isotropic system of macroions and point counterions by Linse,16 where, for monovalent counterions, the finding that Zeff ≈ 10rm/λB is in good agreement with the Alexander et al.12 results. Such a scheme might account for the Zeff/Zm ratio for very large values of the ratio Zm/rm as reported by Roberts and coworkers.10 However, the criteria set forth by the abovementioned charge renormalization theories might not be adequate for macroions with low to intermediate values of Zm/rm. A third approach is based not on specific models and assumed forms of the pair interaction but rather on the pair distribution function. Belloni’s alternative suggestion14 is a renormalization model that resembled the Bjerrum approach by defining a cutoff distance at the inflection point of the pair distribution function. Roberts et al.10 suggested that the cutoff distance as that point for which the macroion-counterion pair distribution function attains the value of unity, that is, when the pair interaction changes from an apparent attraction (viz., negative value) to an apparent repulsion. The effective macroion charge is then the bare macroion charge reduced by the charges of the counterions within that distance. It was recently proposed18,19 that the effective macroion charge can be calculated from the distance at which the macroioncounterion interaction energy attains a balance with the thermal energy kBT. The counterions within this distance are said to be bound and reduce the bare macroion charge to give an effective charge. Recently, one of us (K.S.S.) was involved in Brownian dynamics (BD) simulations19 using this thermal criterion to calculate Zeff/Zm in macroion solutions. The general trend of the results is similar to that seen in the experimental data summarized by Roberts and co-workers.10 Because of technical difficulties, the highest value of Zm that could be treated in the BD simulations was Zm ≈ 350. A purpose of the present work is to try to simulate higher values of Zm using canonical MC simulations in conjunction with PB calculations. Apart from the thermal criterion of Schmitz,18,19 the concentration criterion of Roberts et al.10 will be used to calculate the effective macroion charge for a variety of physical conditions, and hence, a detailed comparative study of the two criteria will be performed. It is emphasized that the thermal and concentration criteria differ from other charge renormalization models by including microions that are found far from the macroion surface. The rest of the paper is organized as follows. In the next section, the system geometry and the details of the MC simulations and PB calculations are discussed. The results are reported in section 3, and in section 4, some general conclusions are drawn. 2. Methods The system under examination is a spherical macroion of radius rm surrounded by its charge cloud under salt-free conditions. Hence, |Zm/Zc| ) Nc is the number of counterions, all with common radius rc. The macroion volume fraction for these calculations is thus
φm )
Vm (4π/3)rm3 ) V V
(3)
where Vm is the volume of the macroion and V is the total volume of the computation box, which contains only one macroion. The standard Wigner-Seitz spherical cell model is employed in both the MC and nonlinear PB calculations, with the ionic species (18) Schmitz, K. S. Langmuir 1999, 15, 4093. (19) Sanghiran, V.; Schmitz, K. S. Langmuir 2000, 16, 7566.
being treated as rigid ions in a dielectric continuum, viz., the primitive model. The macroion is placed in a spherical cell, the size of which is determined by the volume fraction φm of the macroion. The (fixed) macroion charge Zm|e| is spread evenly over the macroion surface leading to a uniform surface charge density of σ ) Zm|e|/(4πrm2). The radius of the spherical cell, rcell, in terms of the macroion radius rm, is given by
rcell )
( ) 1 φm
1/3
rm
(4)
The concentration of the counterions is calculated from the macroion free volume, cc ) Nc/Vf, where the free volume, Vf, is given by
Vf ) V - (4π/3)(rm + rc)3 ) V - (4π/3)rs3
(5)
Here, rs ) rm + rc is the radius of the exclusion sphere centered on the macroion and represents the distance of closest approach of a counterion to the macroion. Thus, the center of a counterion can traverse the region (rm + rc) e r e rcell encompassing the volume Vf. Notice that this definition of the free volume differs slightly from the generally used definition of V(1 - φm). Admittedly, the volume of the counterions is not considered in calculating the free volume, as this correction is negligible compared to the correction due to the macroion volume. For instance, for Nc counterions, we have
4 Nc π rc3 3
(
)
[34 π(r
m
]
Nc rm 1+ rc
(
)
+ rc)3
)
3
which is on the order of 10-4 at the highest value of Nc (1000) for the particle sizes considered in this work. The molar counterion concentration cc is based on the macroionfree volume of the cell and is related to the macroion charge and radius through
cc ) )
|Zp/Zc| NAVf
[
Nc
NAV 1 - φm
(
)]
rc 1+ rm
3
(6)
where NA is Avogadro’s number. Note that, upon taking the limit rc f 0, eq 6 reduces to the form given in ref 10 for one species of counterions in the cell. Defining ψ(r) as the mean electrostatic potential at a distance r from the origin (center of the macroion) and assuming one species of counterions is present, the thermal criterion18,19 would imply that the counterions within a radius r < rtherm contribute to the effective charge where
|Zceψ(r)rtherm)| )1 kBT
(7)
In an analogous manner, the concentration criterion of Roberts et al.10 stipulates that counterions within r < rconc contribute to the effective charge where now the macroion-microion pair distribution function satisfies the condition
g(r)rconc) ) 1
(8)
2.1. Monte Carlo Simulations. The standard Metropolis algorithm was employed in a canonical ensemble MC simulation of the spherical cell model with the model (physical) parameters being defined by eqs 3-6. Although the actual number of configurations sampled varied depending on the macroion charge and volume fraction, in general, the first ∼107 configurations were used for equilibration, while statistics were obtained over the following ∼108 configurations. The displacement parameter
Effective Charge of a Particle in a Salt-Free Suspension (for the random move of a particle) was adjusted until an acceptance rate of around 50% was obtained. Also, in creating a histogram for low macroion charges (Zm e 50) it was seen that a relatively bigger bin size (than that for Zm > 50) led to better statistics and less noise in the MC data. As a test of the MC numerics, we were successful in reproducing some known MC simulation results in the literature for the osmotic pressure of spherical micellar solutions.20,21 The thermal radius, rtherm, and the corresponding Zeff were evaluated in two essentially different ways. In the first, the mean electrostatic potential, ψ(r), was constructed from its formal definition (see, for example, the review in ref 22) in terms of the positions of the individual particles at any given configuration, with the ensemble average of ψ(r) then being obtained by averaging over the configurations as the MC run proceeded. In the second method, Poisson’s equation, viz.
∇2ψ(r) ) -
|e| Z n g (r) 0r c c MC
(9)
where gMC(r) is the MC macroion-microion pair distribution function (obtained at the end of the MC run), was integrated using a fourth-order Runge-Kutta procedure.23 It proved useful here to smooth the MC data by using a least-squares polynomial fit. The differences in the results of the two procedures were very small, usually within statistical errors. The evaluation of rconc and the corresponding Zeff was more straightforward. The function gMC(r) was simply examined to locate the value of r ) rconc where gMC ) 1. Again, the process benefited from the use of a smoothed gMC(r). 2.2. Nonlinear Poisson-Boltzmann Calculations. These calculations involve the numerical solution of the PB equation in a spherical geometry appropriate for the spherical cell model containing one species of counterion only. The equation was used earlier in studies of the structure and thermodynamics of micelles with or without an added salt in the cell model.20,21,24 Monte Carlo simulations20,21,24 have shown that, for monovalent, small, simple microions and/or not-too-dense a solution, the PB approach gives fairly adequate results. As most of the present calculations involve monovalent counterions only at low to moderate concentrations, the PB results are expected to be reasonable. Although the counterions are uncorrelated point ions in this approximation, we assign, in the Stern25 spirit, a distance of closest approach rs (cf. eq 5) of a counterion to the central macroion. Thus, in the PB theory, the macroion-counterion radial distribution function gmc(r) can be written as
gmc(r) ) H(r-rs) gmc(r)rcell) exp{-Zc|e|β[ψ(r) - ψ(r)rcell)]} (10) where H(x) is the Heaviside function. The mean electrostatic potential, ψ(r), is governed by the PB equation, which is simply eq 9 with gMC(r) replaced by the above gmc(r). For convenience only, we put ψ(r)rcell) ) 0. This is a routine practice in the literature (see, for example, ref 21). Using eq 7, the thermal or energy criterion can be expressed in terms of the PB g(r)’s as
loge
[
]
gmc(r)rtherm) gmc(r)rcell)
)1
(11)
Hence, the effective macroion charge Zeff is determined by (20) Wennerstro¨m, H.; Jo¨nsson, B.; Linse, P. J. Chem. Phys. 1982, 76, 4665. (21) Bhuiyan, L. B.; Outhwaite, C. W.; Bratko, D. Chem. Phys. Lett. 1992, 193, 203. (22) Outhwaite, C. W. In Statistical Mechanics (Specialist Periodical Reports); The Chemical Society: London, 1975; Vol. 2, p 188. (23) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. Numerical Recipes, The Art of Scientific Computing; Cambridge University Press: New York, 1986. (24) Amaro, M. L.; Bhuiyan, L. B.; Outhwaite, C. W. Mol. Phys. 1995, 86, 725. (25) Stern, O. Z. Elektrochem. 1924, 30, 508.
Langmuir, Vol. 18, No. 11, 2002 4213 Zeff ) Zp + 4πccZc
∫
rcriterion
rs
gmc(r)r2 dr
(12)
where the upper limit is rcriterion ) rtherm (or rconc).
3. Results and Discussion The MC simulations were performed up to a maximum macroion charge of Zm ) 1000 and for a wide range of physical and model parameters. These parameters were chosen in keeping with earlier works (see, for example, refs 5, 18, and 19). Thus, most of the calculations were done for monovalent counterions and a Bjerrum length of around 7.14 Å with varying macroion volume fraction and species’ radii. We have, however, extended the calculations to macroion systems involving divalent counterions at two different macroion volume fractions. The PB equation was solved by a quasi-linearization iteration technique.26 The method had been used earlier with much success in solving the PB and modified Poisson-Boltzmann equations in the bulk22 and electric double layers including the spherical cell model.21,24 3.1. Macroion-Counterion Pair Distributions. We begin this discussion by studying some macroion-counterion pair distribution functions given in Figure 1 for monovalent counterions at Zm ) 100 (subfigure A), Zm ) 500 (subfigure B), and Zm ) 1000 (subfigure C). Results are presented at φm ) 0.01 (curves a), 0.03 (curves b), and 0.05 (curves c) for each Zm. As can be expected, with decreasing φm, the contact value of gmc(r) increases because of the relatively lower screening of the macroion charge, with the PB predictions being very close to the MC data. Analogous observations were recorded by Stevens et al.5 in their MC work for Zm ) 450, φm ) 0.04 (cf. Figure 2 of ref 5). Such results, in the so-called weak interaction regime, are not surprising; for example, even at the highest Zm and φm values considered in Figure 1, the macroion surface charge density and the counterion concentration are 0.01 C/m2 and 0.0146 M, respectively, which are modest values. It is, of course, well-known in the (bulk) electrolyte and electric double layer literature that the PB theory affords useful descriptions of systems containing (only) monovalent ions at dilute to moderate solution concentrations (e0.5 M) and/or low surface charge densities in the case of a double layer. This is because the interionic correlations, neglected in the mean-field approximation, are quite small for monovalent ions at low solution densities. 3.2. Energy versus Concentration Criteria for Determining Effective Charge. We next turn to assessing the comparative behavior of the effective macroion charge accirding to the energy and concentration criteria. The results for Zeff/Zm, along with those for the corresponding rtherm/rm and rconc/rm values calculated at rm ) 113 Å, rc ) 1 Å, λB ) 7.14 Å, are shown in Figures 2 and 3. Both (Zeff/Zm)therm and (Zeff/Zm)conc are monotonically decreasing with Zm/rm. The rate of decrease also decreases as the bare macroion charge increases, with the Zeff/Zm plots becoming flatter at higher Zm. However, there are important quantitative differences between the two criteria. At low values of the bare charge, the ratio (Zeff/ Zm)therm tends to unity. This was also observed in the BD simulations.19 Indeed, in the limiting case Zmf 0, the thermal criterion predicts (Zeff/Zm)therm f 1. As the bare charge increases at any given macroion volume fraction, the macroion-counterion interaction energy begins to dominate, and more counterions are attracted to near the (26) Bellman, R.; Kalaba, R. Quasilinearization and Nonlinear Boundary Value Problems; Elsevier: New York, 1965.
4214
Langmuir, Vol. 18, No. 11, 2002
Figure 1. Macroion-counterion pair distribution function gmc(r) for a salt-free spherical colloidal suspension containing monovalent counterions at macroion radius rm ) 113 Å, counterion radius rc ) 1 Å, Bjerrum length λB ) 7.14, and macroion charges of (A) Zm) 100, (B) Zm) 500, and (C) Zm) 1000. In each subfigure, results are shown at macroion volume fractions of (a) φm ) 0.01, (b) φm ) 0.03, and (c) φm ) 0.05. The symbols represent MC data, whereas the continuous lines represent the PB results. The vertical lines on curves a indicate the locations of the thermal criterion radius rtherm (solid line) and concentration criterion radius rconc (dashed line). (Such lines are not shown for curves b and c to avoid clutter.)
macroion surface, resulting in an increase of rtherm/rm and a consequent decrease in (Zeff/Zm)therm. However, with this buildup of counterion concentration close to the macroion surface, there is also a stronger repulsion that a distant diffuse-layer counterion must feel in trying to migrate nearer to the macroion. These two competing effects lead to both rtherm/rm and (Zeff/Zm)therm becoming gradually flatter at higher Zm/rm. The physical situation with regard to the concentration criterion is quite different though. Note that, even at very low bare charge, gmc(contact) > 1 as a result of packing effects, so that (Zeff/Zm)conc does not go to unity as Zm f 0. The PB calculations of Roberts et al.10 and the present calculations support this conjecture. Accordingly, one would expect (Zeff/Zm)therm > (Zeff/Zm)conc for low values of the bare charge, which is seen in Figure 2. However, the roles reverse at higher Zm, and the fact that now (Zeff/ Zm)conc > (Zeff/Zm)therm seems to suggest that, in this regime, the energy criterion is more effective than the concentration criterion in “neutralizing” the bare charge. The relative behavior of rtherm/rm versus rconc/rm is interesting. In dilute solutions (low bare charge), gmc(r) also decays slowly, hence, the high values of rconc/rm (cf. Figure 3). As the bare charge increases, however, gmc(r) begins to fall off more steeply, resulting in a decrease in rconc/rm. In
Mukherjee et al.
Figure 2. Effective charge of a macroion in a salt-free colloidal suspension containing monovalent counterions. The quantity plotted is Zeff/Zm as a function of Zm/rm at macroion radius rm ) 113 Å, counterion radius rc ) 1 Å, Bjerrum length λB ) 7.14, and macroion volume fractions of (A) φm ) 0.01, (B) φm ) 0.03, and (C) φm ) 0.05. In each subfigure, results are shown for (a) the concentration criterion and (b) the thermal criterion. The symbols represent MC data, whereas the continuous and dashed lines represent the PB and screened Coulomb results, respectively. The screened Coulomb results are given for the thermal criterion only.
passing, we remark that the PB results in Figures 2 and 3 continue to be in notably good agreement with the simulations, whereas the screened Coulomb results are displaced from the MC data. 3.3. Extension to Divalent Counterions. Figures 4 and 5 show the results of our extension of the calculations to divalent counterion systems. These calculations were again performed with physical parameters of rm ) 113 Å, rc ) 1 Å, and λB ) 7.14 Å and at macroion volume fractions of 0.005 and 0.01. The results for the effective charge are shown in Figure 4, while Figure 5 shows the rtherm/rm and rconc/rm results. It is readily seen that these results for divalent counterions exhibit the same general qualitative trends as their corresponding monovalent counterparts. However, there are significant quantitative differences between the results for mono- and divalent counterion systems, as shown in Table 1. It is noted, for example, that both rtherm/rm and rconc/rm for the divalent ion simulations are larger than the corresponding ratios for the monovalent counterions for the same values of φm and Zm/rm. This might seem paradoxical as the divalent counterions achieve a higher concentration at the surface of the macroion as compared to the monovalent counterions. However, the thermal criterion is not a surface effect involving only macroion-counterion interactions. Indeed, the thermal criterion extends into the bulk medium and represents the electrostatic energy of the counterion in
Effective Charge of a Particle in a Salt-Free Suspension
Figure 3. Thermal criterion and the concentration criterion radii about a macroion in a salt-free colloidal suspension containing monovalent counterions. rtherm/rm (curves a) and rconc/ rm (curves b) are plotted as a function of Zm/rm at macroion radius rm ) 113 Å, counterion radius rc ) 1 Å, Bjerrum length λB ) 7.14, and macroion volume fractions of (1) φm ) 0.01, (2) φm ) 0.03, and (3) φm ) 0.05. The symbols represent MC data, whereas the continuous lines represent the PB results.
Langmuir, Vol. 18, No. 11, 2002 4215
Figure 5. Thermal criterion and the concentration criterion radii about a macroion in a salt-free colloidal suspension containing divalent counterions. rtherm/rm (curves a) and rconc/rm (curves b) are plotted as a function of Zm/rm at macroion radius rm ) 113 Å, counterion radius rc ) 1 Å, Bjerrum length λB ) 7.14, and macroion volume fraction of (1) φm ) 0.005 and (2) φm ) 0.01. The symbols represent MC data, whereas the continuous lines represent the PB results. Table 1. Comparison of the Zeff/Zm, rtherm/rm, and rconc/rm Results at Constant Zm/rm Ratio and Constant Om for Mono- and Divalent Counterion Systemsa Zm/rm (Å-1)
rtherm/rm MC PB
rconc/rm MC PB
(Zeff/Zm)therm MC PB
(Zeff/Zm)conc MC PB
0.443 1.33 2.65 3.98 5.31 7.08 8.85
1.51 2.06 2.18 2.22 2.23 2.25 2.26
1.51 2.06 2.20 2.23 2.26 2.27 2.28
monovalent counterions 2.81 2.80 0.918 0.919 2.38 2.38 0.569 0.572 2.05 2.06 0.327 0.332 1.88 1.90 0.227 0.231 1.78 1.80 0.173 0.176 1.70 1.71 0.131 0.134 1.62 1.65 0.106 0.108
0.675 0.510 0.344 0.259 0.208 0.165 0.139
0.677 0.516 0.349 0.263 0.213 0.170 0.143
0.443 1.33 2.65 3.98 5.31 7.08 8.85
1.90 2.15 2.18 2.19 2.19 2.19 2.20
1.91 2.19 2.26 2.27 2.28 2.29 2.19
divalent counterions 2.57 2.57 0.720 0.726 2.01 2.06 0.317 0.332 1.74 1.80 0.165 0.173 1.61 1.68 0.110 0.120 1.52 1.61 0.083 0.091 1.46 1.54 0.061 0.069 1.40 1.50 0.049 0.055
0.587 0.333 0.196 0.140 0.111 0.085 0.072
0.598 0.349 0.213 0.155 0.124 0.098 0.082
a Calculations performed at r ) 113 Å, r ) 1 Å, φ ) 0.01, and m c m λB ) 7.14 Å.
Figure 4. Effective charge of a macroion in a salt-free colloidal suspension containing divalent counterions. The quantity plotted is Zeff/Zm as a function of Zm/rm at macroion radius rm ) 113 Å, counterion radius rc ) 1 Å, Bjerrum length λB ) 7.14, and macroion volume fractions of (A) φm ) 0.005 and (B) φm ) 0.01. In each subfigure, results are shown for (a) the concentration criterion and (b) the thermal criterion. The symbols represent MC data, whereas the continuous lines represent the PB results.
the mean electrostatic potential field, in which the distribution of charge is an implicit factor. Whereas more divalent counterions accumulate near the surface, their more concentrated divalent charge compared to the more dilute monovalent charge results in an anticorrelation with subsequent divalent charges at distances away from the surface but close enough to distinguish the charge distribution. Such a distribution might be accompanied by a well-defined inflection point in g(r). Such inflection points are a well-known phenomenon with adsorption
isotherms when those of highly interactive bound ligands are compared with weakly interacting (almost independent) ligands. A highly repuslive interaction between bound ligands results in a larger average separation distance between these ligands when compared to the weakly interacting system at the same degree of saturation in the isotherm. In other words, the highly repulsive interaction effectively excludes certain configurations in the performance of averaging over the desired quantities. One might refer to this effect as an inflationary thermal radius. Another important difference between the mono- and divalent results is that the PB predictions for the divalent counterion systems, although qualitatively agreeing with the corresponding MC simulations, show greater deviations from the MC simulations. This is expected and consistent with the known behavior of the PB theory visa`-vis the presence of divalent ions in electrolyte and
4216
Langmuir, Vol. 18, No. 11, 2002
Mukherjee et al.
Table 2. Universality in the (Zeff/Zm) Results for Constant Zm/rm Ratio and Macroion Volume Fractiona (Zeff/Zm)therm φm
Zm
rm (Å)
rc (Å)
0.01 0.01 0.01 0.04 0.04 0.04 0.04
33 66 132 60 225 35 350
11 22 44 40 150 10 100
1 2 4 0.4 1.5 0.1 1
a
Zm/rm
(Å-1)
3.00 3.00 3.00 1.50 1.50 3.50 3.50
(Zeff/Zm)conc
MC
BD
PB
MC
BD
PB
0.283 0.299 0.310 0.500 0.511 0.221 0.255
0.2643b 0.311b 0.03624b 0.1729b
0.315 0.315 0.315 0.513 0.513 0.260 0.260
0.299 0.317 0.328 0.453 0.464 0.247 0.292
0.240b -
0.337 0.337 0.337 0.466 0.466 0.299 0.299
Calculations performed for monovalent counterions at λB ) 7.18 Å. b Reference 15. Table 3. Comparison of Some MC and BD Results for (Zeff/Zm)therm and (Zeff/Zm)conca (Zeff/Zm)therm φm
0.0001 0.001 0.005 0.01 0.01 0.03 0.04 0.04 0.04 0.05
Zm 66 66 110 66 220 200 225 350 100 110
rm (Å) 22 22 28 22 22 113 150 100 10 28
rc (Å-1) 2 2 2 2 2 1 1.5 1 1 2
Zm/rm (Å-1) 3.00 3.00 3.93 3.00 10.00 1.77 1.50 3.50 10.00 3.93
(Zeff/Zm)conc
MC
BD
PB
MC
BD
PB
0.450 0.353 0.238 0.299 0.120 0.450 0.511 0.255 0.0747 0.241
0.4164b
0.465 0.369 0.252 0.315 0.126 0.454 0.513 0.260 0.110 0.254
0.444 0.359 0.263 0.317 0.183 0.430 0.464 0.292 0.0967 0.278
0.260c 0.240b -
0.459 0.375 0.279 0.337 0.197 0.435 0.466 0.299 0.165 0.299
0.233b 0.305b 0.1547b 0.167b 0.304c 0.311b 0.1729b 0.0173b 0.248b
a Corresponding PB results also shown for completeness. Calculations performed for monovalent counterions at λ ) 7.18 Å. b Reference B 18. c V. Sanghiran and K. S. Schmitz, unpublished results.
polyelectrolyte systems, which stems from the neglect of the fluctuation potential (ion-ion correlations) and ionic excluded-volume effects in the classical mean-field approximation. In fact, a nonlinear least-squares fit of the MC simulations for divalent ions requires at least two screened potential forms to fit the entire radial range from the macroion. Such double potential forms were not needed for the monovalent case except for macroions of high charge. In other words, the potential from the surface of the particle to the boundary of the cell cannot be described by a single potential or a single screening parameter as there are near-field and far-field components.27 3.4. Universality of Zeff/Zm and Comparison with BD Calculations. There are, in the literature, at least two types of universality in regard to macroion systems and effective charge. Shklovskii17 defines universality in the context of charge saturation. In the case of strongly charged spheres, the PB equation predicts that the effective charge becomes independent of the bare charge. However, Shklovskii also notes that this universality disappears when correlations are taken into account. The second use of universality refers to the data summarized by Roberts and co-workers,10 where somewhat of a master plot is obtained when Zeff/Zm is plotted as a function of Zm/rm (cf. Figure 3 of ref 10). These data indicate a monotonically decreasing Zeff/Zm with increasing Zm/rm. The universality issue in Zeff/Zm results was also raised by Sanghiran and Schmitz,19 the idea being that, at a given value of φm for the same Zm/rm ratio, the values of Zeff/Zm should be consistent. It is the second definition that we employ in the present study. The Zeff/Zm plots in Figures 2 and 4 can be considered universal. Indeed, these results are qualitatively in agreement with the experimental data summarized by Roberts et al.10 (cf. Figure 3 of ref 10). Similar features have also been observed in the works by Alexander et al.,12 Stevens et al.,5 and Linse,16 for example. To examine this issue more closely, some Zeff/Zm results are collected (27) Mukherjee, A. K.; Bhuiyan, L. B.; Schmitz, K. S., manuscript in preparation.
in Tables 2 and 3, along with the corresponding BD results19 for comparison. A slight increase in the MC Zeff/ Zm with macroion size is noticeable in Table 2 compared to that for the PB approximation. This is again attributable to the ionic correlations and excluded-volume (size) effects. Note that, in the PB approach, the ions are uncorrelated point charges. Correlation effects are also manifested in the choice of the ratio of radii, rm/rc. We can compare our MC simulations for Zm ) 66, rm ) 22 Å, and rc ) 2 Å with those of Lobaskin, Lyubartsev, and Linse9 for Zm ) 60, rm ) 20 Å, and rc ) 4 Å. In both studies, Zm/rm ) 3 Å-1, so one might expect on the basis of universality that the ratio Zeff/Zm should be the same in each case. However, in our study, we find Zeff/Zm ) 0.299 (φm ) 0.01), whereas Lobaskin et al.9 reported this ratio as Zeff/Zm ) 0.355 (φm ) 0.008). Although the macroion volume fractions used are slightly different, a part of the deviation from universality is attributed to the larger size of the counterions employed by Lobaskin, Lyubartsev, and Linse.9 The distance of closest approach to other counterions and the macroions is larger for the larger counterions given that all other parameters are identical. Thus, the universality is expected to be better satisfied for small values of the ratio rc/rm, viz., very large macroions when the counterion radius is held fixed. Further comparisons of the MC, PB, and BD results for Zeff/Zm for the two criteria are shown in Table 3. The PB theory is again seen to reproduce the MC data quite well. It ought to be mentioned that a (relatively) smaller macroion size and a fairly high macroion charge would imply a higher macroion surface charge density and a higher mean counterion concentration, which, in turn, would lead to greater deviation in the mean-field predictions from the simulations. The MC and BD results are consistently qualitatively in agreement overall and, in some instances, are fairly close. Furthermore, the minima observed in both the MC and PB results for the volume fraction dependence of Zeff/Zm at a given Zm, as shown in Figure 6, were also seen earlier in the BD calculations.19,28
Effective Charge of a Particle in a Salt-Free Suspension
Langmuir, Vol. 18, No. 11, 2002 4217 Table 4. Comparison of Present MC Results for (Zeff/Zm) with Those Using the Criteria of Alexander et al.8 and Stevens et al.5 a (Zeff/Zm)therm
(Zeff/Zm)conc
Zm
φm
Zeff/Zmb
Zeff/Zmc
MC
PB
MC
PB
250
0.01 0.05 0.2 0.4
0.332 0.384 0.544 0.720
0.232 0.232 0.248 0.292
0.487 0.600 0.806
0.487 0.485 0.600 0.808
0.460 0.424 0.472 0.501
0.460 0.448 0.476 0.509
500
0.01 0.05 0.2 0.4
0.352 0.214 0.338 0.506
0.246 0.125 0.152 0.180
0.270 0.276 0.358 0.513
0.274 0.280 0.364 0.515
0.298 0.312 0.368 0.423
0.302 0.316 0.374 0.433
a PB results also shown for comparison. Calculations (for monovalent counterions) performed with the parameters used by Stevens et al.,5 rcell ) 150 Å, rc ) 3 Å, and λB ) 7.14 Å. b Criterion of Alexander et al.8 as calculated by Stevens et al.5 c Criterion of Stevens et al.5
Figure 6. Macroion volume fraction dependence of the effective charge of a macroion in a salt-free spherical colloidal suspension containing monovalent counterions obtained at bare macroion charge Zm ) 60, macroion radius rm ) 60 Å, counterion radius rc ) 2 Å, and Bjerrum length λB ) 7.18 using (a) the thermal criterion, and (b) the concentration criterion. The symbols represent MC data, whereas the continuous lines represent the PB results.
It is noted however, from the data in Tables 2 and 3, that the BD results are systematically lower than the MC results, apart from the three exceptions in Table 3 where the results are comparable, with the MC being only slightly lower than the BD values. This systematic difference stems from the differences in the underlying philosophy of the two methods of simulation. First, all of the moves are accepted in the BD simulations, except, of course, when a move might result in an overlap of hard spheres. In contrast, the MC simulations have an energy criterion regarding acceptance if the new energy is above some threshhold value (again the overlap of hard spheres falls into this category). The second difference is that the BD simulations are described as a biased MC simulation.29 The biased nature of the BD simulations means that there is a directional movement dictated by the cumulative force on a particle, viz., the total Coulombic force on it due to the rest of the particles in the system. Thus, there will be a strong tendency of the counterions to remain in the vicinity of the macroion because of this additional force not present in the MC simulations, resulting in the need for long simulation times to reach stability. Closely tied to this is another factor that might have also contributed to the lowering of the BD effective charge, which is the “technical difficulty” mentioned earlier,18 namely, the relatively small time window allowed on the available computer facility at the time of these simulations. In some of these cases, only 106 configurations were used in the averaging process. Indeed, longer simulation times resulted in more stable values.28 3.5. Comparison with Other Criteria for Calculating Effective Charge. The pressure criterion of determining the renormalized charge advocated by Alexander et al.12 in their pioneering work of on the subject (cf. Introduction), later adopted by Stevens et al.5 in their simulations, is compared with the energy and concentration criteria in Table 4. The calculations (for monovalent (28) Sanghiran, V. Ph.D. Thesis, University of Missouri-Kansas City, Kansas City, MO, 2001. (29) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Clarendon Press, Oxford Science Publications: Oxford, U.K., 1993.
counterions) were performed at the same physical parameters used by Stevens et al.,5 viz., rm ) 150 Å, rc ) 3 Å, and λB ) 7.14 Å. It is interesting to notice that the Alexander et al.12 results for Zeff/Zm underestimate the thermal criterion results but overestimate the concentration criterion results at Zm ) 250 (φm ) 0.2, 0.4), but at Zm ) 500, both criteria are underestimated at higher φm. The relative behaviors of these two criteria follow the pattern seen in Figure 2; for example, at the lower Zm/rm ratio, (Zeff/Zm)therm > (Zeff/Zm)conc, but at the higher Zm/rm ratio, (Zeff/Zm)conc > (Zeff/Zm)therm. Our results are thus qualitatively consistent with these results from the literature. The differences in the results in the table can be traced to the differences in the methodologies used in the renormalized charge calculation. For example, in the thermal criterion, the value of the thermal radius first increases and then decreases as the ratio Zm/rm increases. The approach of Alexander et al.12 was to assume that the strong Coulombic interaction of the counterions at the surface of the macroion rendered them bound to the macroion and that this energy was balanced by the entropy of mixing. Stevens et al.5 showed that the effective charge calculated by the method of Alexander et al.12 was consistently larger than that obtained from their simulations. However, in the method of Stevens et al., the pressure was calculated from the macroion-counterion distribution at the farthest “shell” from the macroion. Thus, these two studies represent the extremes in the distribution function in the calculation of the effective charge of the macroion. In the present MC simulations, both the energy and concentration criteria involve the counterion distributions. (The mean electrostatic potential is really a weighted function of the distributions.) We have calculated an effective charge (energy criterion) based on the relative energetics of the counterions. If the total Coulomb interaction experienced by a counterion is larger than or equal to the thermal energy, then that counterion is said to be a contributor the effective charge of the macroion. Similarly, the concentration criterion utilizes some intermediate value of gmc(r) between its values at the macroion surface and at the cell boundary. However, importantly, all of the results in Table 4 have common qualitative features, viz., the decrease of Zeff/Zm with Zm/ rm, and the increase of Zeff/Zm with φm. Groot30 proposed a “condensation” model for spherical particles patterned on the Manning condensation model for linear polyions (see, for example, the review in ref 31) (30) Groot, R. D. J. Chem. Phys. 1991, 95, 9191. (31) Manning, G. S. Q. Rev. Biophys. 1978, 11, 179.
4218
Langmuir, Vol. 18, No. 11, 2002
Mukherjee et al.
and, in the process, obtained an effective charge based on a combination of PB theory and MC simulations. Some comments on the appropriateness of a Manning-type condensation for spherical geometries are in order. As pointed out by Zimm and Manning,32 the physical picture that accompanies the mathematics is the bringing from an infinite distance of discrete charges that eventually comprise the linear array of charges making up the linear polyelectrolyte. More specifically, the electrostatic interaction in the Manning model is limited to the phosphate groups along the DNA backbone. The role of counterions is to neutralize these phosphate-phosphate interactions while maintaining a statistical equilibrium with the bulk solution via the entropy of mixing. Thus, both the summations of the equally spaced phosphate groups and the entropy of mixing terms have logarithmic forms, which when combined lead to the concept of condensation to prevent the argument from going to infinity. Le Bret and Zimm33 have shown that condensation of the Manning type does not occur in a spherical geometry, although it does for planar systems. That is, in the case of spherical macroions, the counterions can be diluted away in the limit of infinite dilution. Whereas the Manning model does not have a term corresponding to the Coulomb interaction of a counterion with the discrete line charges, the spherical model of Groot explicitly states such an interaction. In this regard, the assumption of Groot that Zmax ≈ πrm2/λB2 ) 35 has seemingly no basis in either theory or concept. The effective charge of Groot was obtained from his eq 26,30 viz.
rcell3 - rm3
Zeff ) 4π
rcell3
-b
3
∫br
cell
F(r)r2 dr
(13)
where F(r) is the counterion density profile and the parameter b is chosen on the basis of a uniform distribution of counterions to give the bare macroion charge, which corresponds to the limit λB/rc f 0 and Zm f 0. In other words, the prefactor 4π(rcell3 - rm3)/(rcell3 - bm3) was fixed for all simulations. The parameters used in his simulations were (in units of the particle diameter) rm ) 0.5, b ) 0.65, and rcell ) 0.7. It is clear from the relative values of these numbers that the range of the computational cells available for the counterions, and hence the calculation of the effective charge, is quite small. The observation of a maximum in the value of Zeff might not necessarily be due to the association of counterions with the macroion but rather to the mutual repulsion of the counterions from the thin shell of relative thickness 0.2. In contrast, the values of rtherm and rconc in our simulations varied with the bare charge and volume fraction of the macroion. Thus, for a highly charged macroion, a constant value of the effective charge might be due to the compensating features of an increase in the number of counterions in the vicinity of the macroion and a reduction in the value of rtherm or rconc. 4. Concluding Remarks The main achievement of this paper has been an indepth MC study of two recently proposed criteria for calculating the effective charge of a spherical macroion in a charge-stabilized colloidal suspension. Evidently, from the above discussion, the determination of the effective charge depends on the criterion one brings to bear. The properties of the energy and concentration criteria ex(32) Zimm, B. H.; Manning, G. S. J. Chem. Phys. 1965, 43, 4250. (33) Le Bret, M.; Zimm, B. H. Biopolymers 1984, 23, 287.
amined in this paper have many similarities and are also qualitatively similar to some other such criteria in the literature, for example, in the decrease of Zeff/Zm with Zm/ rm and the appearance of a minimum in the macroion volume fraction dependence of Zeff/Zm. The former feature is also borne out in experimental measurements of effective charge. It is of interest to observe that the concept of thermodynamically bound and free ions has been used in linear polyelectrolyte solutions.34 A significant feature of the present calculations is the universality of the effective charge ratio Zeff/Zm for the same values of Zm/rm for realistic macroion parameters in which the macroion radius is considerably larger than the radius of the microions. This universality is most readily apparent in the PB results in Table 2, where the microions are assumed to be point charges. As indicated by the MC simulations, deviations in the ratio Zeff/Zm are less apparent for the smaller values of the ratio rc/rm. The observation and verification of the notion that the ratio Zeff/Zm is dependent only on the ratio Zm/rm means that one can use simulation techniques to investigate highly charged colloidal systems (Zm g 104) in terms of more manageable equivalent systems. For example, simulations for Zm ) 10 000 and rm ) 500 ought to give the same results for Zeff/Zm as in the case of Zm ) 1000 and rm ) 50. The implication is that the distribution functions have the same shape and characteristics. As regards which method of determining the effective charge is “best” of the ones available, we remark that the concentration criterion underestimates the effective charge at low Zm/rm as it does not predict a limiting value of Zeff/ Zm ) 1 as Zm/rm f 0. Other criteria in the literature have implied constraints in their implementation procedures. Extrapolation to either the surface of the macroion or the boundary of the computation cell has the implicit assumption that the variation in the profile is relatively slow and smooth. Such profiles are expected for low values of the ratio Zm/rm. Furthermore, such extrapolation procedures are to regions of high uncertainty. For example, the concentration of counterions can be very low near the cell boundary, giving rise to large uncertainties in the evaluation of the concentration. In contrast, the Belloni method14 focuses on regions of the distribution function where the uncertainties in gmc(r) can be considered to be unimportant. However this method has the requirement of a sharply changing gmc(r) to obtain a well-defined inflection point and is probably useful for large values of the ratio Zm/rm. The proposed method of the thermal radius criterion has the advantage of the Belloni method in regard to the uncertainties in gmc(r) but does not require an inflection point. Furthermore, it is important to note that the thermal radius criterion is applicable to clusters of macroions as well as isolated macroions. That is, within a thermal radius, two or more macroions can reside. For example, in the case of two charged spheres in close proximity, the hypernetted chain calculations of Sa´nchez-Sa´nchez and Lozada-Cassou35 indicate that the counterion distribution totally engulfs both charged spheres. It was recently reported that BD simulations coupled with the thermal radius criterion likewise indicated that the counterions were to be associated with both charged spheres.36 Hence, in contrast to the other models discussed in this paper, the thermal radius method is sensitive to macroion positions within clusters. Such clusters, with relatively (34) Pack, G. R.; Wong, L.; Lamm, G. Biopolymers 1999, 49, 575. (35) Sa´nchez- Sa´nchez, C.; Lozada-Cassou, M. Chem. Phys. Lett. 1992, 190, 202. (36) Schmitz, K. S. Langmuir 2001, 17, 8028.
Effective Charge of a Particle in a Salt-Free Suspension
high potentials because of the macroions present, should draw in excessive numbers of counterions relative to that predicted by isolated macroion parameters. The thermal radius method should then be sensitive to such effects, for then the number of free counterions, which contribute to the conductivity of the supporting medium, would be reduced when compared with the predicted numbers of free ions. Such effects have not been part of the focus of the present paper, which was designed primarily to examine the concept of the thermal radius and its validity by comparison with other studies on the topic of the effective charge. As indicated previously, our results are in reasonably good agreement with the MC simulations of Lobaskin, Lyubartsev, and Linse.9 Hence, any thermodynamic calculations and applications valid for other models likewise are applicable to the thermal radius approach advocated in our study. The present study therefore adds to the credibility that the thermal radius can provide valuable information about the distribution
Langmuir, Vol. 18, No. 11, 2002 4219
of the microions in such systems and the influence of the microion distribution on thermodynamic and structural properties. To the best of our knowledge, the theoretical and simulation work on effective charge in charge-stabilized colloidal suspensions has, to date, involved monovalent counterions only. Our results for divalent counterions clearly show the relative effects of the increased interionic correlations on the system properties. We hope to consider such cases involving multivalent counterions using formal statistical mechanical approaches and simulations in a future publication. Acknowledgment. Partial support of an NSF Grant 0137273 is acknowledged. We are grateful to Professor C. W. Outhwaite for comments that helped to improve the paper. L.B.B. acknowledges an internal grant through FIPI, University of Puerto Rico. LA011625T
