ARTICLE pubs.acs.org/Langmuir
Interactions Between Charged Surfaces with Ionizable Sites Stephen A. Barr* and Athanassios Z. Panagiotopoulos Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544, United States ABSTRACT: A key factor controlling the interactions between surfaces in aqueous solutions is the surface charge density. Surfaces typically become charged though a titration process where surface groups can become ionized based on their dissociation constant and the pH of the solution. In this work, we use a Monte Carlo method to treat this process in a system with two planar surfaces with explicitly described ionizable sites in a salt solution. We focus on a system with a surface density of ionizable sites set to 4.8 nm2, corresponding to silica. We find that the surface charge density changes as the surfaces come close to contact due to interactions between the ionizable groups on each surface. In addition, we observe an attraction between the surfaces above a threshold surface charge, in good agreement with previous theoretical predictions based on uniformly charged surfaces. However, close to contact we find the force is significantly different than for the uniformly charged case.
’ INTRODUCTION Surfaces are typically composed of ionizable surface groups which become charged through a titration process when in contact with a solvent. This situation is ubiquitous in nature, and understanding the electrostatic interactions between surfaces is therefore of great importance. Over the past decades, charged surfaces have been studied extensively with the electrostatic interactions typically calculated using mean-field theories based on the Poisson Boltzmann (PB) equation, such as DLVO theory.1,2 This works well when the electrostatic interactions are weak and correlations between the ions are not important. Recently, a new theory has been developed, called strong coupling theory (SC), which does take into account counterion correlations and predicts attraction between like charged surfaces, which DLVO theory fails to predict,3,4 in good agreement with simulation results. These theories were originally developed assuming a uniform charge distribution on the surfaces; however, there has been some work in which the discrete nature of charged surfaces is taken into account for both SC5,6 and mean-field theories using the PB equation.7,8 In addition to theoretical approaches, charged surfaces have also been extensively studied using computer simulations with both uniform911 and heterogeneous5,1215 surface charge distributions. One important factor controlling the interactions between charged surfaces is the surface charge density, σ. While many theoretical and simulation studies take σ as an input to the calculations, there has also been extensive work in understanding and being able to predict the surface charge density based on the solution conditions. For this, the surface complexation model has traditionally been used,1618 which is based on the Poisson Boltzmann equation. Since this is a mean-field theory, correlations are neglected. The effect of correlations on ionization has been the subject of much study,19 and recently, a grand canonical titration (GCT) Monte Carlo method has been developed to calculate σ based on the solvent conditions, such as the salt r 2011 American Chemical Society
concentration, the pH, and the intrinsic dissociation constants of the ionizable surface groups.20 Using this method, Labbez et al. have demonstrated the importance of interactions between the ionized surface groups on the resulting value of σ for a single surface.21 In particular, they have compared this method with the surface complexation model and with Monte Carlo simulations of a uniformly charged surface and found that both of these underestimate σ for systems with divalent salt. GCT, which treats the surface groups explicitly, agrees well with experimental results of a silica surface. This method has also been successfully applied to cement pastes22 and clay particles.23 The GCT method has also been used to study a system of two parallel surfaces. In this case, an increase in σ close to contact is observed for large pH values, which is qualitatively different than the results of the surface complexation model which instead predicts a monotonic decrease in σ as the surface separation decreases.21 In this paper, we use this method to gain further insight into the interactions between two charged surfaces close to contact. We study how σ varies as the separation between the surfaces is decreased, focusing on the effect of correlations between the ionizable surface groups and the counterions. We also look at the pressure between the surfaces and compare this with results for a system with a uniform charge density in order to highlight the role of correlations between the ionizable surface groups on opposing walls.
’ MODEL AND METHODS Our system consists of two parallel surfaces which are constructed from discrete ionizable sites. The surfaces are immersed in a solution with divalent salt, CaCl2. We use the slab geometry where the system is Received: April 13, 2011 Revised: June 2, 2011 Published: June 06, 2011 8761
dx.doi.org/10.1021/la201353u | Langmuir 2011, 27, 8761–8766
Langmuir
ARTICLE
periodic in only the x and y directions. A snapshot of the system is shown in Figure 1. All particles in the system, including the ionizable sites on the surfaces, interact through the purely repulsive Weeks-Chandler-Andersen (WCA) potential,24 8 2 ! !6 3 12 > > d 5 > 4 d < þ ε if rij < 2ð1=6Þ d 4ε rij rij VWCA ðrij Þ ¼ ð1Þ > > > ð1=6Þ :0 d if rij > 2 with ε/kB = 298 K and d = 0.4 nm. The electrostatic interaction between two particles with charges, qi and qj, a distance rij apart is given by the Coulomb potential qi qj V ðrij Þ ¼ 4πε0 εr rij
ð2Þ
where ε0 is the permittivity of free space and εr is the relative dielectric constant of the solvent, which is modeled implicitly. These interactions are calculated using the particleparticleparticle-mesh Ewald method of Hockney and Eastwood.25 Because the system is only periodic in two dimensions, we add a vacuum layer in the nonperiodic direction of at least 3 times the system width and add a dipole correction term.2628 In order to model the ionization of the surface sites, we use the method outlined by Labbez and J€onnson.20 In this method, each ionizable site has an intrinsic dissociation constant, K0, and the titration
Figure 1. Snapshot of the system showing two surfaces separated by a salt solution. The top and bottom surfaces are constructed out of ionizable sites on a square lattice with the ionized sites shown in light blue and deionized sites shown in gray. The positive salt ions in solution are shown in red and the negative ions are blue, to match the negatively charged ionized surface sites.
can be described as Si OH h Si O þ H þ ,
K0 ¼
aSiO aH aSiOH
ð3Þ
where ai is the activity of species i. The sites are able to change their charge state through a grand canonical Monte Carlo titration move at constant pH and chemical potential of the salt. Deprotonation can be thought of as a three-step process: the addition of a salt molecule, the deprotonation of the surface site and release of a proton, and the exchange of the proton and the negative ion with the bulk. However, because the counterions in our system are divalent, two randomly selected sites must simultaneously deprotonate to maintain charge neutrality. These steps can be combined into a single move which is accepted or rejected according the standard Metropolis acceptance criterion26,29,30 " accðo f nÞ ¼ min 1,
Np ðNp 1ÞV exp β ΔU μCa ðNCa þ 1ÞðNd þ 1ÞðNd þ 2Þ
þ 2 ln 10 pH pK0
#
ð4Þ
where Np is the number of protonated sites, Nd is the number of deprotonated sites, NCa is the number of Ca ion in the old configuration, μCa is the chemical potential of the Ca ions, V is the volume, and ΔU is the difference in energy between the old configuration and the trial configuration. The expression for protonation is given by " Nd ðNd 1ÞNCa exp β ΔU þ μCa accðo f nÞ ¼ min 1, V ðNp þ 1ÞðNp þ 2Þ # ð5Þ 2 ln10 pH pK0 For a titration move, we choose randomly between protonation and deprotonation with equal probability. The salt is included using a grand canonical Monte Carlo scheme in which the chemical potential is fixed and the number of particles fluctuates. To maintain electroneutrality, we insert or delete neutral groups of ions. Since we are using CaCl2, a group consists of two ions with charge 1e, representing the chlorine, and one ion with charge þ2e to represent the calcium ion. μ is determined for both salt ions independently using the Widom insertion method31 in a bulk salt solution at a fixed concentration. The chemical potential of the Ca ion is also used for grand canonical Monte Carlo titration moves. For a grand canonical move, we choose randomly between insertion and deletion with equal probability. We also choose randomly between grand canonical moves of the salt and titration moves with equal probability and a single run consists of 2 106 Monte Carlo moves for equilibration followed by
Figure 2. Ionization fraction, R, as a function of surface separation, h, for different values of pH pK0. In the left figure, εr = 78.4, and in the right figure, εr = 39.2. Error bars are smaller than the symbol sizes. 8762
dx.doi.org/10.1021/la201353u |Langmuir 2011, 27, 8761–8766
Langmuir
ARTICLE
Figure 3. Snapshot from a typical simulation showing a bridging configuration. The Ca ion, shown in red, is bridging the two ionized sites, shown on the top and bottom in light blue. The deionized surface sites are shown in gray. Figure 5. Fraction of Ca ions that are in a bridging configuration with ionized sites on both walls, fbridging, as a function of pH pK0 for ε r = 78.4 and εr = 39.2. Error bars are smaller than the symbol sizes.
Figure 4. Pair correlation function, g(r), between the ionized surface sites and the Ca ions for εr = 78.4 and pH pK0 = 1. We use the location of the first minimum, r = 0.68 nm, as a cutoff distance when determining the fraction of Ca ions in a bridging configuration. This location is similar for all systems investigated. 2 107 moves. A typical simulation takes 50 h on a single 3.00 GHz Intel Core 2 Duo processor. In all the simulations, the temperature is set to T = 298 K and the ionizable sites have a pK0 of 7.7 to match the pK0 of silica.16,21 These are placed on a square 400 400 lattice with a surface density of 4.8 nm2. For a single surface, this model has been shown to produce surface charge densities in good agreement with experimental results for a silica surface up to pH = 9.21
’ RESULTS A key property which controls the interactions between surfaces is the surface charge density. This will, in general, vary as the separation between the surfaces is decreased, since interactions between ionizable sites on opposing walls become important at small separations. In Figure 2, we show the degree of ionization in a 100 mM CaCl solution for εr = 78.4, corresponding to water, and with εr = 39.2 which corresponds to a polar solvent such as dimethylformamide that enhances the effective strength of the electrostatic interactions. At large separations, the walls do not interact with each other and behave as single isolated surfaces. At separations closer than 1 nm, however, deviations from the bulk behavior begin to occur. With εr = 78.4, there is an increase in ionization when pH pK0 > 1, while at pH pK0 = 1, the degree of ionization decreases close to contact. The increase in ionization for larger values of the pH close to contact occurs because of a bridging configuration, where a divalent counterion is in close proximity to multiple ionized sites on both walls. This
mitigates the mutual repulsion between the ionized sites and is an energetically favorable configuration. Figure 3 shows an example of this bridging. At pH pK0 = 1, however, there are fewer ionized sites, and while the bridging configuration is energetically favorable, it is entropically unfavorable. With only a small fraction of sites ionized, it is unlikely for two sites on opposing walls and a Ca ion to line up, making bridging configurations less likely to occur. This can be made more precise by determining the fraction of Ca ions that are in a bridging configuration, fbridging. A Ca ion is counted as being in a bridging configuration if it less than 0.68 nm from an ionized site on both the top and the bottom surfaces. We chose 0.68 nm as the cutoff for this calculation because it is the location of the first minimum in the pair correlation function between Ca ions and ionized sites, shown in Figure 4. fbridging is shown in Figure 5 for the closest separation as a function of pH pK0 for εr = 78.4 and εr = 39.2. In the case where the ionization fraction decreases at contact, pH pK0 = 1 and εr = 78.4, fbridging is only 0.09. With εr = 39.2 at the same pH, the electrostatic interactions are stronger, and even though the degree of ionization is similar to εr = 78.4, fbridging is almost 5 times higher, and in this case where entropy is less important, there is indeed an increase in ionization fraction close to contact. At larger pH values, the ionization fraction is larger making bridging more likely, and at pH pK0 = 2, nearly all of the Ca ions are in a bridging configuration. We also observe nonmonotonic behavior in σ for some systems at h ≈ 0.6 nm, with a maximum in σ at low H and a minimum at high pH . At this separation, two distinct layers of Ca ions begin to form. In this crossover regime, between two strongly interacting surfaces sharing a monolayer of counterions and well-separated walls, the effect the interacting layers of counterions have on σ is unclear. In addition to the ionization fraction, we also determined the pressure, P, between the surfaces by calculating the z component of the forces acting on the ionizable sites and normalizing by the surface area A ! Nsites 1 Nsites DUi, WCA DUi, Coulomb þ ð6Þ P¼ A i¼1 Dz Dz i¼1
∑
∑
where Nsites is the number of ionizable sites on one wall and Ui, WCA and Ui,Coulomb include the interactions of all particles acting on site i. For the uniformly charged system, the second 8763
dx.doi.org/10.1021/la201353u |Langmuir 2011, 27, 8761–8766
Langmuir
ARTICLE
Figure 6. Pressure as a function of separation, h, for both the ionizable system and uniformly charged system: (top left) pH pK0 = 1; (top right) pH pK0 = 0; (bottom left) pH pK0 = 1; (bottom right) pH pK0 = 2. The insets show the corresponding correlations between ionized sites on opposing walls. Error bars for the pressure are smaller than the symbol sizes.
term in eq 6 is replaced by Nions
qσ
σ
2
i þ ∑ 2ε0 εr i ¼ 1 2ε0 εr
ð7Þ
to account for the electrostatic interactions. In Figure 6, we show the pressure as a function of separation for different values of the pH and compare this with simulations in which the charge is constant and uniformly distributed on the surface. Uniformly charged systems have been shown to give different results for σ than ionizable systems for a single surface.21 Because of this, for the uniformly charged system we set σ to the average surface charge density obtained from simulations of the ionizable system. In this way, we are able to directly compare the effects of siteion and sitesite correlation on the pressure. As in the ionizable system, the salt is included grand canonically. At larger separations, the uniform system and the ionizable system yield identical results for the pressure; however, at closer separations there are deviations. These differences can be explained by looking at correlations between ionized sites on opposing walls, shown in the insets of Figure 6. We measure this correlation by calculating the fraction of bridging configurations, fsim fsim ¼ Nbridging =Nsites
ð8Þ
where Nbridging is the number of occurrences in which two opposing sites are both ionized and Nsites is the number of ionizable sites on one surfaces. This is then normalized by frandom, which is the fraction of bridging configurations which would occur if the same number of ionized sites were randomly distributed on each surface Nionized 2 ð9Þ frandom ¼ Nsites
where Nionized is the number ionized sites on one surface. frandom is then equal to the probability that a site on one surface is ionized, Nionized/Nsites, and the probability that the opposing site is also ionized, assuming an equal number of ionized sites on both surfaces. Because the number of ionized sites varies throughout the simulation, the quantity fsim/frandom is computed for each frame and then averaged over the entire simulation. We see in the insets of Figure 6 that at close separations there are positive correlations in all cases, meaning there are more instances of two corresponding sites on opposing walls both ionized than if they were randomly placed. These result from a bridging configuration and will act to reduce the pressure, and indeed, we see in these cases that the pressure is lower for the ionizable system than for the uniform case for all pH values. Even though, for larger values of the pH, the value of the positive correlation is small at the closest separation, we see from Figure 5 that nearly 100% of these have a counterion bridging them, which leads to the large difference in pressure between the ionizable and the uniform cases. At slightly larger separations, for pH pK0 > 0 anticorrelations occur, and in these cases, the pressure of the uniformly charged system is more attractive than the ionizable system. At larger separations, these correlations decay and the uniformly charged system has the same pressure as the ionizable system. Because ions are included explicitly in both the ionizable and the uniform simulations, the differences are only due to siteion and sitesite correlations, as ionion correlations are present in both systems. We also see a small difference in the position of the attractive minimum between the ionizable and the uniformly charged systems for pH pK0 > 0. While the range of the attraction is expected to be determined by the ion size,32 because we are using explicit surface sites which have a soft interaction with the ions as opposed to smooth walls with hard interactions, the extent to 8764
dx.doi.org/10.1021/la201353u |Langmuir 2011, 27, 8761–8766
Langmuir
ARTICLE
ions is included.3335 Because we only have purely repulsive interactions between counterions, we do not observe this phenomenon; however, the interplay between raft formation and ionization would be a rich problem for future study. It would also be interesting to extend this model to study systems such as spherical particles to investigate the effects of ionizable sites on colloidal stability.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
Figure 7. Two-dimensional pair correlation function, g2d(r), for a range of surface separations, h, for the ionizable system with pH pK0 = 1. The inset shows g2d(r) for the corresponding uniformly charged system.
which this shift is a real effect is difficult to determine within the resolution of our data. Lateral correlations between the counterions are also considered. In Figure 7, we compare the two-dimensional pair correlation function, g2d(r), in the plane parallel to the walls for a range of surfaces separations in the ionizable systems with pH pK0 = 1 and the corresponding uniformly charged system. For the closest separation, we consider all counterions, and for larger separations, we only consider counterions within the layer closest to the surface within a width of 0.35 nm, equal to the width of the closest surface separation. We see long-range order for the closest separation. Increasing the surface separation to 0.55 nm, however, results in much less long-range order, and the structure is the same as wellseparated surfaces. We also observe that the uniformly charged system is nearly identical to the ionizable system, which indicates that correlations between the ionized sites do not significantly impact the counterion structure.
’ CONCLUSIONS In this work, we have used a grand canonical Monte Carlo method that accurately reproduces surface charge densities to gain insights into the interactions between charged surfaces. We found that the surface charge density varies with the surface separation and can either increase or decrease depending on the interplay between energy and entropy, energy favoring an increase in ionization due to favorable bridging interactions between the counterions and the ionized sites on opposing walls. These bridging configurations result from the surface groups responding to their local environment. This also affects the pressure, and we find that a uniformly charged system overestimates the pressure because that model does not take into account sitesite or siteion correlations. The uniform approximation is good, however, at larger separations where these correlations are not important. The ionizable method, on the other hand, is able to accurately determine both the surface charge density and the interaction between the surfaces, even in cases where correlations are important. In this work, we focused on a system with a relatively high surface density of ionizable sites. In this case, the charged state of a site is strongly influenced by its neighboring sites, and lowering the surface density could lead to qualitatively different behavior, which could be studied in future work. In similar systems, ionic raft formation has been observed when an attraction between the
’ ACKNOWLEDGMENT This work was supported by grant DE-SC0002128 from the Department of Energy, Office of Basic Energy Sciences. ’ REFERENCES (1) Derjaguin, B. V.; Landau, L. Acta Physicochim. USSR 1941, 14, 633–662. (2) Verwey, E. J.; Overbeek, J. T. G. Theory of the Stability of Lyophobic Colloids; Elsevier: Amsterdam, 1948. (3) Naji, A.; Netz, R. R. Eur. Phys. J. E 2004, 13, 43–59. (4) Naji, A.; Jungblut, S.; Moreira, A. G.; Netz, R. R. Physica A 2005, 352, 131–170. (5) Moreira, A.; Netz, R. Europhys. Lett. 2002, 57, 911–917. (6) Jho, Y.; Park, G.; Chang, C.; Pincus, P.; Kim, M. Phys. Rev. E 2006, 73, 021502. (7) Lukatsky, D.; Safran, S.; Lau, A.; Pincus, P. Europhys. Lett. 2002, 58, 785–791. (8) Henle, M.; Santangelo, C.; Patel, D.; Pincus, P. Europhys. Lett. 2004, 66, 284–290. (9) Guldbrand, L.; J€onsson, B.; Wennerstrom, H.; Linse, P. J. Chem. Phys. 1984, 80, 2221–2228. (10) Moreira, A. G.; Netz, R. R. Eur. Phys. J. E 2002, 8, 33–58. (11) Pellenq, R.; Caillol, J.; Delville, A. J. Phys. Chem. B 1997, 101, 8584–8594. (12) Khan, M.; Petris, S.; Chan, D. J. Chem. Phys. 2005, 122, 104705. (13) van megen, W.; Snook, I. J. Chem. Phys. 1980, 73, 4656–4662. (14) Akesson, T.; Jonnson, B. J. Chem. Phys. 1985, 89, 2401–2405. (15) Messina, R.; Holm, C.; Kremer, K. Eur. Phys. J. E 2001, 4, 363–370. (16) Hiemstra, T.; van Riemsdijk, W.; Bolt, G. J. Colloid Interface Sci. 1989, 133, 91–104. (17) Sverjensky, D. Nature 1993, 364, 776–780. (18) Borkovec, M. Langmuir 1997, 13, 2608–2613. (19) Travesset, A.; Vangaveti, S. J. Chem. Phys. 2009, 131, 185102and references therein. (20) Labbez, C.; J€onsson, B. Lect. Notes Comput. Sci. 2007, 4699, 66–72. (21) Labbez, C.; J€onsson, B.; Skarba, M.; Borkovec, M. Langmuir 2009, 25, 7209–7213. (22) Pochard, I.; Labbez, C.; Nonat, A.; Vija, H.; J€onsson, B. Cem. Concr. Res. 2010, 40, 1488–1494. (23) Delhorme, M.; Labbez, C.; Caillet, C.; Thomas, F. Langmuir 2010, 26, 9240–9249. (24) Weeks, J. D.; Chandler, D.; Andersen, H. C. J. Chem. Phys. 1971, 54, 5237–5247. (25) Hockney, R. W.; Eastwood, J. W. Computer Simulation Using Particles; IOP Publishing: Bristol, 1988. (26) Frenkel, D.; Smit, B. Understanding Molecular Simulation, 2nd ed.; Academic: San Diego, 2002. (27) Yeh, I.; Berkowitz, M. J. Chem. Phys. 1999, 111, 3155–3162. (28) Yeh, I.; Berkowitz, M. J. Chem. Phys. 1999, 110, 7935–7942. 8765
dx.doi.org/10.1021/la201353u |Langmuir 2011, 27, 8761–8766
Langmuir
ARTICLE
(29) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087–1092. (30) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids; Clarendon: Oxford, 1987. (31) Widom, B. J. Chem. Phys. 1963, 39, 2808–2812. (32) Solis, F. J.; Olvera de la Cruz, M. Phys. Rev. E 1999, 60, 4496–4499. (33) Loverde, S. M.; Olvera de la Cruz, M. J. Chem. Phys. 2007, 127, 164707. (34) Loverde, S. M.; Solis, F. J.; Olvera de la Cruz, M. Phys. Rev. Lett. 2007, 98, 237802. (35) Velichko, Y. S.; Olvera de la Cruz, M. J. Chem. Phys. 2006, 124, 214705.
8766
dx.doi.org/10.1021/la201353u |Langmuir 2011, 27, 8761–8766
