Molecular Solvent Model of Spherical Electric Double Layers: A

Sep 24, 2009 - Molecular Solvent Model of Spherical Electric Double Layers: A Systematic Study by Monte Carlo Simulations and Density Functional Theor...
0 downloads 10 Views 431KB Size
Download PDF
13980

J. Phys. Chem. B 2009, 113, 13980–13987

Molecular Solvent Model of Spherical Electric Double Layers: A Systematic Study by Monte Carlo Simulations and Density Functional Theory Chandra N. Patra* Theoretical Chemistry Section, Chemistry Group, Bhabha Atomic Research Centre, Mumbai 400 085, India ReceiVed: August 12, 2009; ReVised Manuscript ReceiVed: September 05, 2009

The structure of spherical electric double layers is studied using Monte Carlo simulation and density functional theory by considering solvent as the third component. In this molecular solvent model (MSM), ions and solvent molecules are considered as charged and neutral hard spheres, respectively, having equal diameter. The macroion is considered as an isolated hard sphere having uniform surface charge density surrounded by the electrolyte and the solvent. The theory is partially perturbative as the hard-sphere contribution to the one particle correlation function is evaluated using suitably averaged weighted density, and the ionic part is obtained through a second-order functional Taylor expansion around the bulk electrolyte. The Monte Carlo simulations have been performed in a canonical ensemble. The system is studied at varying concentrations of electrolytes, and the solvent molecules, at different valences of the electrolyte, at different macroion radii, and at varying surface charge densities. The theory is found to be in good agreement with the simulation results over a wide range of parametric conditions. The excluded volume effects due to the molecular nature of the solvent are shown to have much richer features in diffuse layer phenomena like layering and charge inversion. Introduction Understandingphenomenaorprocessesoccurringatsolid-liquid interfaces, involving especially charged particles near charged surfaces have been the subject of a wide variety of research in view of its applications in biological, electrochemical, and colloidal systems.1 The stability of colloidal solutions and their coagulation properties are known to depend largely upon various macroscopic conditions, which has direct relevance in many application areas such as dispersions of globular proteins, dendrimers, polyelectrolytes, aerosols, etc.2 Also, the interfacial behavior of polyelectrolyte systems are directly influenced by the presence of salt or the distribution of dissolved electrolyte as can be envisaged in biological macromolecules like DNA, RNA, proteins, etc.3 and in many technologically important industries like paint and emulsions.4 There are various experimental techniques that probe into the interfacial region directly or indirectly like small angle neutron scattering, small-angle X-ray scattering, ellipsometry, surface force measurements, electrophoresis, titration, etc. to give the structural information up to a molecular level precision.5,6 Apart from these experimental techniques, there have been an upsurge in theoretical and computational studies of interfacial phenomena.7 The interfacial region between a charged surface along with the neutralizing diffuse layer of counterions is often referred to as the electric double layer (EDL),8,9 which plays the major role in different aspects of the interfacial phenomena. EDLs, although extensively studied in different geometries,10-17 are relatively less explored in spherical geometries.17,18 The theoretical description of EDL is mostly based on models. In the most preliminary model, called the restricted primitive model (RPM),8 the electrolytic components are considered as charged hard spheres of equal diameter and the solvent is treated as the dielectric continuum. The RPM never considers the solvent as a component, thus the first level of complexity has been brought in the molecular solvent model (MSM),19 where solvent * Electronic mail: [email protected].

molecules are introduced as neutral hard spheres. This model in its minimal form has been quite successful in describing the oscillatory nature of the solvent and the ions around a charged surface, including charge inversion phenomena at a lower overall concentration. Further insight into details of EDL properties with the molecular nature of solvent, such as the polarizing properties, can be gained by considering more “civilized” models, e.g., solvent as a hard sphere with embedded point dipoles20,21 or quadrupole tensors, etc.22 The first and conventional analytical theory used to describe EDL in a spherical geometry is the classical Poisson-Boltzmann (PB) equation, which completely neglects the very important interionic steric correlations.23 Although the theoretical predictions of PB were aceptable for simple systems of low-charged colloids in highly dilute univalent electrolyte solutions,24 for highly coupled systems, the theory can never make predictions regarding the nonmonotonous ion distribution that may lead to phenomena like layering and overcharging.25 Hence, to incorporate the interionic correlations, various liquid state theories like the modified Poisson-Boltzmann approach (MPB),26 integral equation methods (IET),27,28 density functional theory (DFT), etc.17,18 have emerged as more accurate alternatives. The density functional approach has been found to be a highly powerful tool to study29 the structure of electric double layers (EDLs) and, in many instances, has been found to be more accurate than other theories. The theory has been successfully applied to EDLs in different geometries, viz., planar,30-32 cylindrical,14,33,16,34 and spherical double layers.17,18 In a conventional DFT, the grand potential functional is expressed35 as a functional of density distribution that attains a minimum value at equilibrium. The exact functional form of the grand potential functional, in general, is not known for nonuniform systems; however, for certain systems, it is known for a corresponding uniform counterpart. On the basis of the known functional form of uniform systems, a number of different approaches can be developed, perturbatively or nonperturbatively, to obtain the density distributions and free energy expressions for nonuniform

10.1021/jp907790t CCC: $40.75  2009 American Chemical Society Published on Web 09/24/2009

MSM of Spherical Electric Double Layers

J. Phys. Chem. B, Vol. 113, No. 42, 2009 13981

fluids. DFT has been widely applied to study the structure of various nonuniform fluids such as simple liquids,36,37 liquids under confinement,38,39 ionic and dipolar fluids,40,41 and polymers and polyelectrolytes.42,43 Simulations give exact solutions for a given model and, thus, can be used to test the effectiveness as well as the applicability of specific versions of any theory.44 Extensive Monte Carlo simulation data are available for planar10 and cylindrical EDLs;11,45 however, the data of spherical EDLs is scant,46-48 mainly because of the difficulty in incorporation of long-range interaction potentials to account for the boundary effects. An elaborate study of spherical EDLs in a fully asymmetric situation was carried out in recent years.49 However, the incorporation of the molecular nature of the solvent through MSM would be a phenomenally computer-expensive task. There are few reports where Monte Carlo simulations have been performed using the molecular nature of the solvent in planar50-52 and cylindrical EDLs53 emphasizing the role of solvent-induced ordering of electrolyte ions. However, any systematic study of spherical double layers under MSM has not been attempted so far. In the present work, we carry out a detailed investigation of the molecular solvent model (MSM) of spherical double layers through Monte Carlo simulation and density functional theory. The hard sphere contribution is approximated through a weighted density approach, whereas the residual Coulombic interaction is calculated perturbatively around a uniform fluid. The second-order direct correlation function required as input is taken from the mean spherical approximation (MSA). The simulations have been performed in a canonical ensemble where a finite size of the solvent molecules is included. The theoretical predictions of the ionic and solvent density profiles as well as mean electrostatic potentials are found to be in good agreement with the Monte Carlo simulation data under various macroscopic parameters such as salt and solvent concentrations, ionic valences, macroion radii, and charge densities on the macroparticle surface. The importance of hard-core exclusion effects due to the incorporation of neutral hard spheres of solvent on the oscillations in density profiles and magnitude of potential profiles have been reflected in all systems studied here. It is shown that phenomena like charge reversal become stronger due to the influence of steric effects imposed by molecular solvents. Theoretical Formulation Molecular Model. We consider the model for a spherical double layer which comprises an isolated charged spherical macroparticle surrounded by an MSM electrolyte solution. The macroparticle is a charged hard sphere of radius R having a uniform surface charge density Q given as

Q)

Ze 4πR2

(1)

where zR is the valence of ion R, ε is the dielectric constant of the medium taken as 78.358, corresponding to water at temperature, T ) 298 K, and r is the interionic distance. In essence, the model considered for the present study is more suitable for aqueous systems with a cosolvent represented by hard spheres. The ions as well as the solvent interact with other solvent spheres only through hard-sphere repulsive interactions. Similarly, the macroion-ion potential is defined as

{

4πR2QzRe , r > R + σ/2 UMR(r) ) εr ∞, otherwise

where r is the radial distance of the ion from the macroion surface Here again, the interaction of solvent with the wall is purely repulsive of the type of hard-wall hard-sphere interactions. Unless otherwise stated, the diameter of the small ions (σ) is taken as 0.425 nm, the typical value used in most EDL simulations. Density Functional Theory. The classical density functional theory35 starts with an expression for the grand free energy, Ω, as a functional of the singlet density profiles, FR(r), of each of the species. At equilibrium, the grand free energy is minimal with respect to variations in the density profiles and this condition is used to determine the density profiles and the free energy. Without delving into much detail, we present here the relevant equations for the nonuniform density distribution of small ions along the radial direction r, thus obtained as

FR(r) ) F0R exp{-β0zRψ(r) + c(1)hs R (r;[{FR}]) 0 (1)el (1)el 0 c(1)hs R ([{FR}]) + cR (r;[{FR}]) - cR ([{FR}])} (4)

for r > R + σ/2 and zero otherwise. Here, ψ(r) is the mean electrostatic potential due to the macroion surface charge density as well as internal ionic distributions and is given by

ψ(r) )

4πe ε

(

2

∫r∞ ∑ zRFR(r') r' - r'r R

)

dr'

{

zRzβe2 , r>σ uRβ(r) ) εr ∞, otherwise

(2)

(5)

(1)el and c(1)hs R (r;[{FR}]) and cR (r;[{FR}]) represent, respectively, the hard sphere and the electrical contribution to the first-order correlation function. In the present work, we employed the weighted density approximation (WDA) based on the work of Denton and Ashcroft (DA)54 to calculate c(1)hs R (r;[{FR}]), by approximating the same with that of uniform fluid counterpart, c˜, at a weighted density, Fj(r), i.e.

c(1)hs ˜ (1)hs(F¯ (r)) R (r;[{FR}]) ) c Within the MSM, the electrolyte solution consists of charged and neutral hard spheres representing small ions and solvent molecules, respectively, having equal diameter, σR () σ), with R representing the components (R ) 1 and 2 represent counterions and coions and R ) 3 represents solvent). The ion-ion pair interaction potential is given by

(3)

(6)

where

F¯ (r) )

∫ dr'hs(|r - r'|;F¯ (r))[ ∑ FR(r')]

(7)

R

and whs(r) is the weight function obtained through the DA prescription.54 However, the electrical contribution has been obtained as perturbation around a uniform fluid given as

13982

J. Phys. Chem. B, Vol. 113, No. 42, 2009

Patra

2

(1)el 0 c(1)el R (r;[{FR}]) - cR ([{FR}]) )

∑ c˜Rβ(2)el(|r -

β)1

r'|;{F0R})(Fβ(r') - Fβ0 )

(8)

(2)el with c˜Rβ taken from the analytical expression within the mean spherical approximation (MSA).55 Thus the density distributions of ions and the solvent molecules are obtained by solving, iteratively, eq 4, and the mean electrostatic potential is calculated using eq 5. Monte Carlo Simulations. Canonical Monte Carlo (CMC) simulations (N,V,T) have been performed using standard Metropolis sampling procedure56 for spherical double layers using the solvent as the third component. The simulation cell is a cubic box with the macroion fixed at its center, surrounded by the small ions and solvent molecules, with the minimum image convention for the ions and solvent molecules imposed on the simulation box, whose length L was considered sufficiently large to neglect macroion-macroion interactions and to obtain stable density profiles. The ions and solvent molecules are considered as charged and neutral hard spheres of equal diameter (σR ) σ). The number of ions for each species was adjusted to satisfy the electroneutrality condition over all the system: N1Z1e + N2Z2e + Ze ) 0. The hard spheres are inserted randomly in the simulation cell and moved through a small distance randomly. The selection criteria for any move is decided through change in total potential energy of the system that comprises macroion-ion, UMR, and ion-ion Coulombic interactions, URβ. We carried out a few simulations for concentrated systems with Ewald sumation for Coulombic interactions. However, we are unable to see any significant change due to Ewald contribution in such a large simulation box. This merely corroborates the fact already ascribed by Degre`ve et al.57 and Degre`ve and Lozada-Cassou,58 that it is unnecessary to carry out Ewald treatment of the Coulombic interaction for spherical EDL systems owning a plasma parameter59-61 Γ < 10 when a sufficiently large simulation box is used. Precisely, in our study, it never crosses the prescribed Γ value. The final equilibrated density distributions are obtained after 8 × 108 total moves, and a final average is done over 5 blocks each having 8 × 108 moves. The acceptance ratio is kept between 0.2 and 0.4 because of the presence of excess steric crowding by highly dense solvent hard spheres.

Results and Discussion This work significantly differs from all the earlier works17,18,62 as it incorporates the molecular nature of the solvent as an individual component in predicting the density profiles of ions and the mean electrostatic potential profiles (MEP) in spherical double layers (SDL). The density profiles of ions and solvent molecules and the mean electrostatic potential profiles of SDL around a spherical macroion of R ) 1.5 nm and Q ) 0.102 cm-2 (corresponding to macroion charge Z/e ) 18) containing 1:1 electrolyte has been compared for RPM and MSM in Figures 1 and 2, respectively. The bulk electrolyte is kept at 1 M concentration in both the models, whereas in MSM, the bulk density of the solvent is taken as F30* ) 0.4. It is amply clear that there is a significant amount of oscillations in the density profiles of counterions and coions in the presence of molecular solvent (MSM) as compared to the case where the solvent is considered as a continuum (RPM). While, monotonic decay of counterion profile and a corresponding growth of coions from the surface is observed for RPM, a distinct peak is observed in

Figure 1. Density profiles of various components for 1 M 1:1 salt around a spherical macroion of R ) 1.5 nm and Q ) 0.102 cm-2 for RPM (0) and MSM (O). Symbols are simulation results, and lines represent theoretical predictions (solid line for MSM, dashed line for RPM). Empty symbols describe the counterion profiles, filled symbols correspond to the coion profiles, and symbols with crosses in them show solvent profiles. The inset shows the expanded view to emphasize the layered structure in MSM.

Figure 2. Mean electrostatic potential profiles for 1 M 1:1 salt around a spherical macroion of R ) 1.5 nm and Q ) 0.102 cm-2 for RPM (0) and MSM (O). Symbols are simulation results, and lines represent theoretical predictions (solid line for MSM, dashed line for RPM).

counterion, coion, and solvent density profiles at a distance of σ from the macroion surface in MSM. The inset in Figure 1 shows the enlarged view of the peak in the density profiles of counterion, coion, and solvent in MSM. There is also an overall increase in the density of both counterions and coions at the macroion surface in the case of MSM over those of RPM. This layering and oscillation properties in density profiles can be better understood in terms of interplay of electrostatic interactions between the macroion and the small ions, the interionic correlations among the small ions, and the hard-core exclusion by the macroion, the small ions, and the solvent molecules. In RPM, the electrostatic interactions with the surface charge leads to accumulation of counterions and depletion of coions from the macroion surface. When solvent is introduced (MSM) into the system, since the density of solvent molecules is much more than the small ions (∼9 times higher when F30*) 0.4 and F10*) F20* ) 0.046 235 9 at 1 M ionic concentration), the hard-sphere repulsions between the ions and the solvent become predominant due to efficient packing leading to preferential accumulation of spheres at the surface. Since this packing contribution is in addition to the electrostatic interactions with the surface charge, this will lead to higher accumulation of counterions and a smaller increase in coion density at the surface. Thus, MSM shows an overall higher concentration of counterions and coions at the surface than RPM. The mean electrostatic potential profiles (MEP) (in Figure 2) show that the magnitude of average potential due to the surface charge decreases in MSM as compared to the RPM. This decrease in magnitude of the potential is merely a consequence of a larger accumulation of counterions at the

MSM of Spherical Electric Double Layers

Figure 3. Density profiles of various components for 1 M 1:1 salt around a spherical macroion of R ) 1.5 nm and Q ) 0.102 cm-2 for MSM having solvent density, F0* 3 ) 0.2 (0, dashed line) and 0.6 (4, solid line). Symbols are simulation results, and lines represent theoretical predictions. Empty symbols describe the counterion profiles, filled symbols correspond to the coion profiles, and symbols with cross in them show solvent profiles. The inset shows the expanded view to emphasize the layered structure in MSM.

Figure 4. Mean electrostatic potential profiles for 1 M 1:1 salt around a spherical macroion of R ) 1.5 nm and Q ) 0.102 cm-2 for MSM with F0* 3 ) 0.2 (0), 0.4 (b), and 0.6 (4). Symbols are simulation results, and lines represent theoretical predictions.

surface in MSM leading to more effective screening of the surface as compared to the RPM. The coion concentration being extremely low, they will have very little contribution on the mean electrostatic potential. In both Figures 1 and 2, the agreement between DFT and simulation results has been very good. However, this agreement is much better for RPM than it is for MSM which may be due to a lack of MSA to accurately represent the coupling between hard-sphere and Coulombic interactions, and this deficiency is increased when the contribution of hard-sphere correlations increases phenomenally, i.e. in MSM. The implications of hard-sphere exclusion effects is more clearly visible when the bulk density of the solvent is varied. Figure 3 depicts the density profiles of ions and solvent molecules at solvent densities, F30*) 0.2 and 0.6 keeping the electrolyte concentrations the same at 1 M. The profiles show sharp structures with strong oscillations at high solvent density F30*) 0.6 and little structures at low solvent density F30*) 0.2. Apart from the prominent second layer of density profiles at ∼5σ in high solvent density case (0.6), there are third and fourth layers of accumulation of ions and solvent at ∼6σ and 7σ distances (inset of Figure 3 shows the expanded view of the peak region of the profiles). These third and fourth layers are absent even in the F0* 3 ) 0.4 case (cf. Figure 1). The concentration of counterions in the first layer (at the macroion surface) is also substantially higher at F0* 3 ) 0.6 as compared to the lower solvent density case (0.2). Higher concentrations of counterions in the first layer also effectively screen the surface charge more, thereby reducing the mean electrostatic potential at the surface and leading to its faster damping at higher solvent density (0.6).

J. Phys. Chem. B, Vol. 113, No. 42, 2009 13983

Figure 5. Density profiles of various components for MSM having 1:1 salt around a spherical macroion of R ) 1.5 nm and Q ) 0.102 cm-2 at 0.1 M (0, dashed line) and 2 M (4, solid line) bulk electrolyte concentrations, respectively. The inset shows the expanded view of 2 M salt to emphasize the charge inversion and third layer formation. The key is the same as in Figure 1.

Figure 6. Mean electrostatic potential profiles for MSM having 1:1 salt around a spherical macroion of R ) 1.5 nm and Q ) 0.102 cm-2 at 0.1 (0), 1 (O), and 2 M (4) salt concentrations, respectively. Symbols are simulation results and lines represent theoretical predictions.

Figure 4 shows the MEP profiles of SDL having 1 M electrolyte solution with solvent density F30*) 0.2, 0.4, and 0.6 indicating a decrease in MEP profile with increase in solvent density. Decreasing the solvent density also leads to a more diffuse and thicker double layer. The theory can reproduce the density as well as the MEP profiles quite accurately as predicted by the simulation results. The variation in ionic concentrations is depicted in Figure 5, which shows the density profiles of ions and the solvent at 0.1 and 2 M bulk concentrations of 1:1 electrolyte keeping the solvent density fixed at F30*) 0.4. Strong oscillations in density profiles of counterions, coions, and solvent molecules are observed in 2 M electrolyte system. A small hump is seen in the counterion profile of 0.1 M concentration at a ∼5σ distance which is due to exclusion effects. The first minimum (at ∼4.5σ) becomes deeper, and the first maximum (at ∼5σ) becomes more prominent as the concentration of the electrolyte increases from 0.1 to 2 M. The concentration of coions also increases at the surface when the bulk electrolyte concentration increases (i.e., for 2 M case) which is compensated for by a slight decrease in layering of counterions at the surface. The inset shows the expanded view of the density profiles of the ions and the solvent for 2 M electrolyte system. The counterion and coion density profiles cross each other at the first maximum, and there is a slight increase in coion density over the counterion density beyond first maximum. This observation has significant impact in electrochemical phenomena as it indicates a possibility to observe charge inversion even for a 1:1 electrolyte provided the solvent and electrolyte concentrations are high. It should be noted that under the same parametric conditions, the RPM does not show any charge inversion phenomena, suggesting that the exclusion effects due to solvent molecules play the crucial role in inducing this phenomena. In 2 M electrolyte solution,

13984

J. Phys. Chem. B, Vol. 113, No. 42, 2009

Patra

Figure 7. Density profiles of various components for MSM around a spherical macroion of R ) 1.5 nm and Q ) 0.102 cm-2 having salts of different valences at 1 M bulk concentrations: (a) 1:1, (b) 2:1, (c) 1:2, and (d) 2:2. The key is the same as in Figure 1.

there exists a small hump in the density profiles at ∼6σ indicating the presence of third layer of ions and solvent accumulation. The magnitude of mean electrostatic potential at the surface decreases with an increase in ionic concentrations (due to larger accumulation of counterions at the surface) leading to better screening of surface potential. Figure 6 shows the MEP profiles for 0.1, 1, and 2 M concentrations of 1:1 electrolyte at F30*) 0.4. The plot shows a lower value of MEP at the surface as well as a faster decaying profile with an increase in bulk electrolyte concentration. The double layer becomes thicker and diffuse as the ionic concentration decreases. Whereas for 0.1 and 1 M concentrations, the MEP profile decays to zero monotonically, at 2 M concentration the profile becomes negative before coming back to zero. This is a direct consequence of charge inversion where the mean electrostatic potential changes its sign from attractive to repulsive and then neutral within a short distance from the interface. The DFT predictions and MC results for density and MEP profiles are in close agreement to each other for all concentrations studied here. The diffuse layer characteristics in an EDL are also known to depend largely on the valence of the salt where increasing the valence of the counterions leads to narrowing of the diffuse layer width owing to better screening of the surface potential (as seen in our earlier study18 of SDL with RPM). In Figure 7, the density distributions of all components are shown as a function of distance from the surface of the spherical macroion of R ) 1.5 nm and Q ) 0.102 cm-2 for the SDL containing salts of different valences (i.e., 1:1, 2:1, 1:2, and 2:2). All the cases show strong peaks at ∼5σ due to hard sphere exclusion effects of the solvent. In the case of 1:1, the three density profiles (ions and solvent) merge at ∼6σ, whereas for 2:1 case, the coion and counterion profiles cross each other after ∼5σ with a higher peak of coion than both solvent and the counterions. This phenomena of charge inversion is even more strong in 1:2 salt where the maxima of the coion becomes higher and that of the counterion becomes lower. On comparing the density profiles of 1:2 and 2:2 salt systems, it is observed that the 2:2 system shows stronger charge inversion than 1:2. The density of the divalent coion (2:2) is significantly low at the surface as compared to the monovalent coion (1:2) which is compensated for in increased peak intensity of the coion profile in the second

Figure 8. Comparison of mean electrostatic potential profiles for MSM around a spherical macroion of R ) 1.5 nm and Q ) 0.102 cm-2 having salts of different valences at 1 M bulk concentrations: 1:1 (O, solid line), 2:1 (], dashed line), 1:2 (0, dashed dotted line), and 2:2 (4, dashed double dotted line).

layer in 2:2 as compared to 1:2 salt. Figure 8 shows the mean electrostatic potential profiles of diffuse layer having ions of different valent salts, i.e. 1:1, 2:1, 1:2, and 2:2. As expected, the MEP profile of 1:1 salt decreases monotonically to zero whereas for 2:1, it passes through negative region before going to zero. In case of the 1:2 salt, the MEP profile starts with the negative value at the surface and after decaying increases to reach a zero value at ∼5σ and again takes positive value before coming back to the zero line. The behavior of the 2:2 salt is similar. The distinctly similar and more or less overlapping profiles of 1:2 and 2:2 salts show that the contribution of coions in calculating mean electrostatic potential is negligibly small. For 1:2 and 2:2 salts, the potential becomes negative at the surface due to large excess accumulation of counterions. This oscillation in the potential profile indicates that the effective charge on the double layer keeps changing in different layers due to excess of counterions, followed by coions, and then again counterions in first and second layers in the diffused region, respectively. It should be pointed out that these long lasting oscillations in MEP profiles is seen only in the case of multivalent counterion containing salts (1:2 and 2:2) at high concentrations (1 M) having an excess of solvent (F30*) 0.4). The theoretical predictions are compared with the simulation results for density and MEP profiles. The theory captures all important phenomena of double layers like layering and

MSM of Spherical Electric Double Layers

J. Phys. Chem. B, Vol. 113, No. 42, 2009 13985

Figure 9. Density profiles of various components for MSM for 1 M 1:1 salt around a spherical macroion of R ) 1.5 nm and at varying surface charge density (Q). The symbols (O, solid line), (0, dashed line), and (4, dashed double dotted line) represent values of Q ) 0.102, 0.204, and 0.306 cm-2, respectively. The key is the same as in Figure 1.

Figure 11. Density profiles of various components for MSM for 1 M 1:1 salt around spherical macroions with surface charge density Q ) 0.102 cm-2 and various radii (R). Symbols 4, 0, and O represent R ) 0.5, 1, and 1.5 nm, respectively. The key is the same as in Figure 1

Figure 10. Mean electrostatic potential profiles for MSM for 1 M 1:1 salt around a spherical macroion of R ) 1.5 nm and with surface charge density, Q ) 0.102 (O, solid line), 0.204 (0, dashed line), and 0.306 cm-2 (4, dotted line). Symbols are simulation results, and lines represent theoretical predictions.

oscillations, charge inversion, etc., but the quantitative agreement between the two deviates for higher valent counterion containing salts such as 2:1 and 2:2, especially at the peak region of the density maximum. The theory underestimates the mean electrostatic values especially for higher valent salts. This indeed is a general trend for all theories where the input correlation function is from the MSA. A little improvement in the approach may make the theory more accurate. The effect of variation in the surface charge density on diffuse layer characteristics has been studied by varying Q for a 1 M 1:1 electrolyte solution having solvent density, F0* 3 ) 0.4. Figure 9 shows the density profiles of counterions, coions, and solvent spheres at three different surface charge densities, viz. at Q ) 0.102, 0.204, and 0.306 cm-2 (corresponding to Z/e ) 18, 36, and 54). It is observed that the value of first minimum and first maximum of the counterion layer increases as Q increases. All the density profiles merge together beyond the second layer of accumulation, irrespective of the surface charge. This happens as the hard sphere exclusion effects are similar due to fixed number of spheres, thus oscillations in density profiles are damped within shorter distance. The increased surface charge density is neutralized by increased accumulation of counterions and corresponding depletion of coions in first and second layer of ions at the surface. The mean electrostatic potential profiles for the corresponding systems are shown in Figure 10. The value of MEP at the surface is maximum for Q ) 0.306 and it decreases with decrease in value of Q. The profile decays to zero much faster for high Q values because larger accumulation of counterions neutralize the surface charge within first two layers of accumulation. The diffused layer width increases with decrease in the value of Q indicating a more diffuse and thicker

Figure 12. Mean electrostatic potential profiles for MSM for 1 M 1:1 salt around spherical macroions with surface charge density Q ) 0.102 cm-2 and various radii, R ) 0.5 (4), 1 (0), and 1.5 nm (O). Symbols are simulation results, and lines represent theoretical predictions. The inset shows the profiles from the closest approach to the macroion surface.

double layer. Even though the charge inversion is not seen in 1 M 1:1 electrolyte at such a high Q ) 0.306, the low MEP value and faster decaying plot indicate that further increasing the Q value may cause charge reversal in the diffused double layer. The theory is found to be in overall agreement with the simulation results. Finally, the influence of the macroparticle radius on the counterion layering is shown in Figure 11, where the surface charge density of the spherical macroion is 0.102 cm-2. The height of the first and second layer peaks of counterion increases as the macroion radius increases. The opposite is true for the coion peaks. In the case of the smallest macroion (R ) 0.5 nm), since the accessibility of the surface for all the spheres is quite high, the counterions, coions, and the solvent segregates quite considerably at the surface. However, for the largest macroion (R ) 1.5 nm), the counterions will dominate at the surface because of larger electrostatic attraction. This is also corroborated by the MEP profiles in Figure 12, which indicates a larger magnitude of the surface potential with increase in the radius of the macroion. The inset of Figure 12 clearly shows that the damping of the potential profiles occurs at the same distance from the surface indicating that the thickness of the

13986

J. Phys. Chem. B, Vol. 113, No. 42, 2009

Patra

TABLE 1: Zeta Potentials ζ ) ψ(R + σ/2) in the MSM Model MC

DFT

MC

conc (M)

Q (cm- 2)

0.1 1.0 1.0 1.0 2.0

0.102 0.102 0.204 0.306 0.102

1.0

0.102

-0.02

1:2 electrolyte -0.07

1.0

0.102

0.14

2:1 electrolyte 0.11

1.0

0.102

-0.04

2:2 electrolyte -0.09

R ) 1.5 nm

1:1 electrolyte 2.074 0.556 0.959 1.201 0.181

2.063 0.553 0.939 1.154 0.201

double layer remains constant. The simulation results for the density and the MEP profiles are quite well-reproduced by DFT. One of the most important quantities in the theory of EDL is the zeta potential ζ, which is the value of the mean electrostatic potential at the closest approach to the charged surface, i.e., ζ ) ψ(R + σ/2). This gives a measure of the thickness of the EDL region in which the total charge is locally nonzero. It can be calculated from eq 5 as

ζ)

4πe ε

(

2

∞ r dr ∑ zRFR(r) r ∫R+σ/2 R + σ/2 R

)

(9)

where FR(r) is the ionic density distribution. The zeta potentials for the system of macroion immersed in different electrolyte solutions at three values of the macroion radius and surface charge density and at different concentrations are summarized in Table 1. It can be easily seen that the present results of MC and DFT go hand in hand with each other. This also indicates that, for multivalent counterions, the zeta potentials become significantly lower, yielding even negative values at certain instances. Concluding Remarks A systematic study of structure of the three component model of the spherical electric double layers is carried out using density functional theory and Monte Carlo simulations where solvent is considered as the third component. The model for spherical double layers consist of a charged macroparticle in the presence of charged and neutral hard spheres of ions and the solvent. The model is more realistic than the restricted primitive model because of the incorporation of the molecular nature of the solvent component. The density functional approach is partially perturbative where hard sphere contribution is calculated through weighted density approximations and the residual Coulombic contribution is treated as a perturbation around a uniform fluid with uniform second-order correlation function used from mean spherical approximation. The system is studied by varying the number of parameters such as electrolyte and solvent concentrations, ionic valences, and the size and surface charge density of the macroion. The density profiles and the mean electrostatic potential profiles as obtained through theoretical predictions are compared with the Monte Carlo simulation results and found to be in overall agreement in most cases. The implications of hard-sphere exclusion effects due to the molecular nature of the solvent are evident in oscillations in density profiles, the layering, and the charge inversion phenomena of electric double

0.461

R ) 1 nm

DFT

MC

DFT R ) 0.5 nm

0.482

0.341

0.342

layer. The qualitative trend clearly resembles that of the planar50,63-65 and cylindrical double layers.53 The packing entropic effect of hard spheres of ions and the solvent molecules causes an excessive accumulation of all the components at the surface, thus the oscillations in the density profiles of counterions, coions, and the solvent molecules increase in presence of molecular solvent as compared to the continuum solvent RPM case. On increasing the solvent concentration, the layering and oscillations in density profiles become quite pronounced because of steric interactions. The effect of increasing the ionic concentrations is similar in both the restricted primitive model (RPM) and molecular solvent model (MSM), where the ionic interactions predominate in the growth of counterions in the first layer. The steric interactions also play the crucial role in inducing the charge inversion seen in density profiles of 2 M 1:1 electrolyte where the coion concentration becomes higher than the counterion concentration beyond the first layer. The multivalent counterions interact more with the macroion and thus screen the surface more effectively causing the diffuse layer width to narrow. The charge inversion becomes stronger in the case of multivalent salts (2:2 salts show stronger charge inversion than 2:1). The oscillations and layering of the ions also increase with increase in the valence of the counterions. The higher surface charge density on the macroion leads to higher accumulation of counterions and depletion of coions at the surface. The value of first maximum and first minimum of the counterion layer increases as surface charge density increases although the width of the diffuse layer width remains constant. The first and second layer peaks of counterion increases and that of coion decreases, as the size of the macroion increases. The mean electrostatic potential (MEP) profiles also corroborates the effects due to hard sphere exclusion on the structure of SDL. The magnitude of the MEP value decreases as a whole in the MSM as compared to the RPM due to higher accumulation of counterions at the surface causing better screening of the charge on the macroion surface. An increase in solvent concentration leads to a drop in the MEP value at the surface and a faster decaying MEP profile. Higher ionic concentrations lead to more layering of counterions, thus again causing lowering of the MEP profiles. AT 2 M ionic concentrations, the charge reversal phenomena is visible in the MEP profile where the profile changes its sign from positive to negative before going to zero. Ionic valency plays a significant role in the layering, oscillations, and the charge inversion phenomena in spherical double layers in RPM as well as in MSM. In the present system, 1:2 and 2:2 salts show stronger charge inversion and faster decaying MEP profiles due to better

MSM of Spherical Electric Double Layers screening of surface charges by multivalent counterions, thus leading to narrower diffuse layers. The value of MEP increases and the width of the diffuse layer decreases with increase in the surface charge of the macroion. The increase in size of the macroion also reduces the width of the diffuse layer. The agreement between DFT and MC simulations is found to be quite good except for a few cases where coupling between hard sphere and Coulombic interactions seems to be quite strong, such as systems with high ionic concentrations, with high valence ions, and high surface charge densities. The DFT can predict the oscillations and trends in the density profiles quite satisfactorily, but in some cases, it proves to be inaccurate in predicting the exact maximum and minimum values at the peak region of the density profiles. The disagreement in the density profiles and mean electrostatic potentials near the particle surface may be further improved by using a more advanced formalism for the hard sphere repulsion. The deficiency in theory may arise due to insufficient treatment of residual Coulombic interactions which is calculated as a perturbation over the uniform system. The mean spherical approximation used as input for this perturbative term can be inaccurate for the given system. Although the model captures most of the consequences of the hard-core exclusion effects due to the molecular nature of the solvent, exact treatment of solvent effects may be incorporated by introducing the polar nature of the solvent molecules. Because the electrostatic interaction will be dramatically increased if the dielectric constant is set as one, it will be interesting to see how the DFT works in the “nonprimitive” limit. Also the variation of zeta potentials due to change in different parameters will be another important area of study. At present, work along these directions is in progress in our laboratory and will be reported in future publications. Acknowledgment. The author gratefully acknowledges Swapan K. Ghosh for useful discussions. It is a pleasure to thank Tulsi Mukherjee for his kind interest and constant encouragement. References and Notes (1) Interfacial Processes and Molecular Aggregation of Surfactants; Narayanan, R., Ed.; Springer: Berlin, 2008. (2) Highlights in Colloid Science; Platikanov, D., Exerowa, D., Eds.; Wiley-VCH: Weinheim, 2008. (3) DNA Interactions with Polymers and Surfactant; Dias, R., Lindman, B., Eds; John Wiley: Hoboken, NJ, 2008. (4) Hara, M. Polyelectrolytes: Science and Technology; Dekker: New York, 1993. (5) Israelachvili, J. Intermolecular and Surface Forces; Academic Press: New York, 1992. (6) Imaging of Surfaces and Interfaces: Lipkowski, J., Ross, P., Eds.; Wiley-VCH: New York, 1999. (7) Levin, Y. Rep. Prog. Phys. 2002, 65, 1577. (8) Carnie, S. L.; Torrie, G. M. AdV. Chem. Phys. 1984, 56, 141. (9) Hansen, J. P.; Lo¨wen, H. Annu. ReV. Phys. Chem. 2000, 51, 209. (10) (a) Torrie, G. M.; Valleau, J. P. J. Chem. Phys. 1980, 73, 5807. (b) Torrie, G. M.; Valleau, J. P. J. Phys. Chem. 1982, 86, 3251. (11) Mills, P.; Anderson, C. F.; Record, M. T., Jr. J. Phys. Chem. 1985, 89, 3984. (12) Gonza´lez-Tovar, E.; Lozada-Cassou, M.; Henderson, D. J. Chem. Phys. 1985, 83, 361. (13) Mier-y-Teran, L.; Suh, S. H.; White, H. S.; Davis, H. T. J. Chem. Phys. 1990, 92, 5087. (14) Patra, C. N.; Yethiraj, A. J. Phys. Chem. B 1999, 103, 6080. (15) Patra, C. N.; Ghosh, S. K. J. Chem. Phys. 2002, 117, 8938. (16) Wang, K.; Yu, Y.-X.; Gao, G.-H. Phys. ReV. E 2004, 70, 011912. (17) Yu, Y.-X.; Wu, J.; Gao, G.-H. J. Chem. Phys. 2004, 120, 7223. (18) Goel, T.; Patra, C. N. J. Chem. Phys. 2007, 127, 034502.

J. Phys. Chem. B, Vol. 113, No. 42, 2009 13987 (19) Tang, Z.; Scriven, L. E.; Davis, H. T. J. Chem. Phys. 1992, 97, 494. (20) Carnie, S. L.; Chan, D. Y. C. J. Chem. Phys. 1980, 73, 2949. (21) Blum, L.; Henderson, D. J. Chem. Phys. 1981, 74, 1902. (22) Wei, D.; Torrie, G. M.; Patey, G. N. J. Chem. Phys. 1993, 99, 3990. (23) Tanford, C. Physical Chemistry of Macromolecules; Wiley: New York, 1961. (24) Petris, S. N.; Chan, D. Y. C.; Linse, P. J. Chem. Phys. 2003, 118, 5248. (25) Levin, Y. Annu. ReV. Phys. Chem. 1996, 51, 209. (26) (a) Outhwaite, C. W.; Bhuiyan, L. B. Electrochim. Acta 1991, 36, 1747. (b) Outwaite, C. W.; Bhuiyan, L. B. Mol. Phys. 1991, 74, 367. (27) Gonza´lez-Tovar, E.; Lozada-Cassou, M. J. Phys. Chem. 1989, 93, 3761. (28) Guerrero-Garcı´a, G. I.; Gonza´lez-Tovar, E.; Lozada-Cassou, M.; Guevara-Rodrı´guez, F. de J. J. Chem. Phys. 2005, 123, 034703. (29) Fundamentals of Inhomogeneous Fluids; Henderson, D. Ed.; Dekker: New York, 1992. (30) Patra, C. N.; Ghosh, S. K. Phys. ReV. E 1993, 47, 4088. (31) Boda, D.; Fawcett, W. R.; Henderson, D.; Sokolowski, S. J. Chem. Phys. 2002, 116, 7170. (32) Pizio, O.; Patrykiejew, A.; Sokolowski, S. J. Chem. Phys. 2004, 121, 11957. (33) Patra, C. N.; Yethiraj, A. Biophys. J. 2000, 78, 699. (34) Wang, K.; Yu, Y.-X.; Gao, G.-H.; Luo, G.-S. J. Chem. Phys. 2005, 123, 234904. (35) Evans, R. In Fundamentals of Inhomogeneous Fluids; Henderson, D., Ed.; Dekker: New York, 1992. (36) Tarazona, P. Phys. ReV. A 1985, 31, 2672. (37) Denton, A. R.; Ashcroft, N. W. Phys. ReV. A 1989, 39, 426. (38) Goel, T.; Patra, C. N.; Ghosh, S. K.; Mukherjee, T. J. Chem. Phys. 2004, 121, 4865. (39) Goel, T.; Patra, C. N.; Ghosh, S. K.; Mukherjee, T. J. Chem. Phys. 2005, 122, 214910. (40) Tang, Z.; Scriven, L. E.; Davis, H. T. J. Chem. Phys. 1992, 96, 4639. (41) Patra, C. N.; Ghosh, S. K. J. Chem. Phys. 1997, 106, 2752. (42) Patra, C. N. J. Chem. Phys. 2007, 126, 074905. (43) Patra, C. N.; Chang, R.; Yethiraj, A. J. Phys. Chem. B 2004, 108, 9126. (44) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1987. (45) Goel, T.; Patra, C. N.; Ghosh, S. K.; Mukherjee, T J. Chem. Phys. 2008, 129, 154906. (46) Degre`ve, L.; Lozada-Cassou, M.; Sa´nchez, E.; Gonza´lez-Tovar, E. J. Chem. Phys. 1993, 98, 8905. (47) Degre`ve, L.; Lozada-Cassou, M. Mol. Phys. 1995, 86, 759. (48) Terao, T.; Nakayama, T. Phys. ReV. E 2001, 63, 041401. (49) Guerrero-Garcia, G. I.; Gonzalez-Tovar, E.; Lozada-Cassou, M.; Guevara-Rodriguez, F. de J. J. Chem. Phys. 2005, 123, 034703. (50) Zhang, L.; Davis, H. T.; White, H. S. J. Chem. Phys. 1993, 98, 5793. (51) Boda, D.; Henderson, D.; Patrykiejew, A.; Sokolowski, S. J. Chem. Phys. 2000, 113, 802. (52) Li, Z.; Wu, J. Phys. ReV. E 2004, 70, 031109. (53) Goel, T.; Patra, C. N.; Ghosh, S. K.; Mukherjee, T. J. Chem. Phys. 2008, 129, 154707. (54) Denton, A. R.; Ashcroft, N. W. Phys. ReV. A 1991, 44, 8242. (55) Waisman, E.; Lebowitz, J. L. J. Chem. Phys. 1972, 56, 3086. 1972, 56, 3093. (56) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087. (57) Degre`ve, L.; Lozada-Cassou, M.; Sa`nchez, E.; Gonza`lez-Tovar, E. J. Chem. Phys. 1993, 98, 8905. (58) Degre, L.; Lozada-Cassou, M. Mol. Phys. 1995, 86, 759. (59) Brush, S. G.; Sahlin, H. L.; Teller, E. J. Chem. Phys. 1966, 45, 2102. (60) Card, D. N.; Valleau, J. P. J. Chem. Phys. 1970, 52, 6232. (61) Valleau, J. P.; Cohen, L. K. J. Chem. Phys. 1980, 72, 5935. (62) Bhuiyan, L. B.; Outhwaite, C. W. Condens. Matter Phys. 2005, 8, 287. (63) Patra, C. N.; Ghosh, S. K. Phys. ReV. E 1993, 48, 1154. (64) Tang, Z.; Scriven, L. E.; Davis, H. T. J. Chem. Phys. 1994, 94, 4527. (65) Boda, D.; Chan, K.-Y.; Henderson, D. J. Chem. Phys. 1998, 109, 7362.

JP907790T