Charge regulation of colloidal particles in aqueous solutions
CCharge regulation of colloidal particles in aqueous solutions
Amin Bakhshandeh, a) Derek Frydel, b) and Yan Levin c) Instituto de F´ısica, Universidade Federal do Rio Grande do Sul, Caixa Postal 15051, CEP 91501-970, Porto Alegre,RS, Brazil Department of Chemistry, Federico Santa Maria Technical University, Campus San Joaquin,7820275, Santiago,Chile
We study charge regulation of colloidal particles inside aqueous electrolyte solutions. To stabilize colloidalsuspension against precipitation, colloidal particles are synthesized with either acidic or basic groups on theirsurface. In contact with water these surface groups undergo proton transfer reaction, resulting in colloidalsurface charge. The charge is determined by the condition of local chemical equilibrium between hydroniumions inside the solution and at the colloidal surface. We use a model of Baxter sticky spheres to explicitlycalculate the equilibrium dissociation constants and to construct a theory which is able to quantitativelypredict the effective charge of colloidal particles with either acidic or basic surface groups. The predictions ofthe theory for the model are found to be in excellent agreement with the results of Monte Carlo simulations.The theory is further extended to treat colloidal particles with a mixture of both acidic and basic surfacegroups.
I. INTRODUCTION
Aqueous solutions are of great importance in biologyand chemistry . In many cases such solutions are ionic.The long-range Coulomb interaction between chargedparticles is the mains source of difficulty for exploring thethermodynamics of such systems . Because of organicfunctional groups many surfaces and membranes acquiresurface charge when placed in water. From their interac-tion with water and acid or base these functional groupscan either lose or gain a proton, becoming charged .The amount of charge gained in this process depends onthe pH of solution and the process is known as chargeregulation (CR) . CR is of great importance in col-loidal science, biology, and chemistry , and is respon-sible for the stability of many different systems .The concept of charge regulation was first described byLinderstrøm-Lang and later developed by many other re-searchers . The first quantitative implementationof charge regulation was done by Ninham and Parsegian(NP) who combined the idea of the local chemical equi-librium with the Poisson-Boltzmann theory introducedby Gouy and Chapman sixty years earlier . The fun-damental assumption of the NP theory is that the bulkassociation constants can be used to study proton trans-fer reactions with the surface adsorption sites. Withinthe NP approach the bulk concentration of hydroniumions is replaced by the local density determined self con-sistently by the Boltzmann distribution, c surf H + = c bulk H O + exp( − βφ ) , (1)where β = 1 /k B T and φ the electrostatic surface poten-tial. a) Electronic mail: [email protected] b) Electronic mail: [email protected] c) Electronic mail: [email protected]
The NP theory relies on Poisson-Boltzmann equationwith the CR implemented as a new boundary condition.The two parameters that determine the boundary condi-tion are the equilibrium constant of the chemical reactiontaking place at the surface, and the surface density of thechemical groups. Within the NP model, the surface ishomogeneous, therefore, the model ignores the discretestructure of surface chemical groups. Another assump-tion is that the equilibrium constant is defined in termsof the concentrations of the reacting species rather thantheir activities. The validity of this assumption needs tobe tested, since the concentration of hydronium ions canbe quite large near a charged surface. Finally, the valueof the equilibrium constant at the surface is assumed tobe the same as for the reaction in the bulk. This is clearlyfar from obvious.In this paper we will focus on spherical colloidal par-ticles with acidic and basic surface groups. If a colloidalparticle has N site basic functional groups on its surfacethen the effective surface charge within the NP theory isfound to be σ = K Bulk N site q c a e − βφ π ( a + r ion ) (1 + K Bulk c a e − βφ ) , (2)where a is the colloidal radius and c a is the bulk concen-tration of strong acid. On the other hand if the surfacehas N site acidic groups, the effective surface charge is σ = − N site q π ( a + r ion ) + K Bulk N site q c a e − βφ π ( a + r ion ) (1 + K Bulk c a e − βφ ) , (3)where K Bulk is the bulk equilibrium association constant,which is the inverse of the acid dissociation constant K a , and q is the elementary proton charge. The elec-trostatic surface potential φ must be calculated self-consistently by combining these expressions with themean-field Poisson-Boltzmann equation.The fundamental ingredient of the NP theory is theequilibrium constant. In the original approach the equi-librium constant for the active sites on the colloidal sur-face was assumed to be the same as for the bulk solution, a r X i v : . [ c ond - m a t . s o f t ] A ug however, in the latter works the equilibrium constant wastreated as a fitting parameter. Clearly this is not verysatisfactory, since it does not allow us to explicitly probethe validity of the theory. Although, the NP theory wasa pioneering first step in understanding charge regula-tion in colloidal systems, in the absence of an explicitmodel on which the theory could be tested, the validityof the underlying approximations of the theory remainsunclear. A different approach was recently advocated byBakhshandeh et al. in which a specific model of of associ-ation was used to calculate exactly the bulk equilibriumconstant for acid . The same acidic groups where thenplaced on top of a spherical colloidal particle and thedensity profiles for hydronium cations and correspond-ing anions were calculated exactly – within this model –using Monte Carlo simulations. Knowledge of the exactequilibrium constant allowed us to explicitly compare theresults of simulations with the NP theory. It was foundthat NP approach deviated significantly from the pre-dictions of simulations. For the specific case of acidicsurface groups Ref. then introduced an alternative ap-proach which was found to be in excellent agreement withthe Monte Carlo simulations. The objective of this pa-per is to extend the results of to colloidal particles withbasic surface groups, as well to the particles containinga mixtures of basic and acidic surface groups.We should note that the present theory applies directlyonly to the specific model of chemical association de-scribed below. There are two levels of approximation thatwe use: 1 – the microscopic model of acid-base associa-tion in terms of the Baxter sticky spheres, and 2 – the ap-proximations used to theoretically solve the model. Theadvantage of this two step approach is that the theorycan be tested against an “exact” solution of the micro-scopic model obtained using the computer simulations.This allows us to separate the possible shortfalls of thetheory from those of the microscopic model. If the theoryagrees with the “exact” solution of the model, any short-falls can then be attributed to the microscopic model ofassociation and not to the approximations which had tobe made to solve the model. The disadvantage of suchapproach is that the theory that we develop applies onlyto the specific microscopic model of acid/base equilib-rium and is not generally universal. This, however, isthe problem with any microscopic theory which does notexplicitly take into account all the quantum effects asso-ciated with the charge transfer at the interface. In theabsence of such “complete” theory, we expect that the ap-proach advocated in the present paper will help to shedinteresting new light on the mechanisms of charge regu-lation of nanoparticles and colloidal suspensions, and inparticular on applicability of mean-field theories to studythis intrinsically strong-coupling problem.There are several possibilities for colloidal surface toacquire charge. The acidic functional groups, such ascarboxyl COOH, can become dissociated due to the fol-lowing reactionHA + H O (cid:29) H O + + A – , (4) resulting in a negatively charged surface. Alternatively,basic functional groups, which originally are not charged,can gain protons from hydronium ions and acquire a pos-itive charge, B + H O + (cid:29) H O + HB + . (5)One example of such functional group is amine NH .The paper is organized as follows. In section II, we in-troduce a model of a colloidal particles with sticky sitesand present the details of Monte Carlo simulations. Insection III we show how the equilibrium association con-stant can be calculated for sticky ions. In section IV wepresent a model for a uniformly sticky colloidal particle.In Section V this model is extended to account for dis-crete basic surface groups, and in Section VI and VIIto discrete acidic groups. In Section VIII we considerparticles with a mixture of both basic and acidic surfacegroups and in Section IX we present our conclusions. II. THEORETICAL BACKGROUND AND MONTECARLO SIMULATIONS DETAILSA. Theoretical background
To study charge regulation of a colloidal surface weuse a model of Baxter sticky spheres . The stickypotential was previously used to study gelation in globu-lar proteins , chemical association in weak acid-basereactions , and sticky-charged wall model . In ourmodel, sticky interactions take on a physical interpre-tation of a chemical bond between proton and acid/basegroups. This is not the first time that a sticky interactionis used to model a chemical bond, to capture some aspectof quantum mechanics in an otherwise classical descrip-tion. The idea has been around for some time and reachesback to 1980 in particular, the work of Blum and Her-rera , and Werthaim for directional chemical bond-ing. Sticky interactions continue to this day being animportant part of soft-matter modeling . In the presentwork sticky interactions, and their quantum-chemical in-terpretation, will be used to study charge regulation ofnanoparticles with surface acid and base groups. Theresults obtained, therefore, are only valid within the spe-cific microscopic model. The model, of course, can be ex-tended and modified to represent a different charge regu-lated system. For example, directional chemical bondingcould be introduced by making a sphere sticky in limitedregions. The size and shape of an absorbing moleculecould be changed. These alternatives are not explored inthe present work.The hydronium ion can become adsorbed on an acidicor basic functional group to form a molecule H + A – or H + B, respectively. To model the binding between H + and A – or B we use an attractive square well potentialwith a repulsive hard core . The first component ofthe interaction potential is the hard-core repulsion, u hs ( r ) = (cid:40) ∞ , r < d, , r > d, (6)where d is the diameter of particles. The second compo-nent is a narrow attractive well, u well ( r ) = , r < d, − ε, d < r < d + ∆ , , r > d + ∆ . (7)To generalize the model, we also include a soft poten-tial u sf ( r > d ), so that the total pair potential becomes u tot = u hs + u well + u sf , u tot ( r ) = ∞ , r < d, − ε + u sf ( r ) , d < r < d + ∆ ,u sf ( r ) , r > d + ∆ . (8)We note that the soft potential u sf ( r ) is effective from r = d , as illustrated in Fig. (1). For illustration, we r u (r) well potentialsoft potentialhard-core r u (r) FIG. 1: Schematic representation of thesticky-hard-sphere potential plus soft interaction. Thesticky part is represented as a narrow well potential; (a)depicts hard-core plus attractive well and soft potentialseparately, and (b) as a combination. An importantobservation is that the soft interaction is active withinthe attractive well.consider a very simple scenario comprised of two particlesinteracting via a pair potential in Eq. (8) and confinedto a spherical region of radius R . To make the demon-stration even simpler, one particle is fixed at the originand only a second particle is free. The resulting partitionfunction has two parts, Z = 4 π (cid:90) Rd r e − β [ u well ( r )+ u sf ( r )] dr =4 πe βε (cid:90) d +∆ d r e − βu sf ( r ) dr + 4 π (cid:90) Rd +∆ r e − βu sf ( r ) dr. (9) Assuming a small ∆, we can expand Z in ∆, yielding Z = 4 πd e βε e − βu sf ( d ) (cid:20) ∆ (cid:18) − e − βε (cid:19) + (cid:18) d − dβu sf ( r ) dr (cid:12)(cid:12)(cid:12)(cid:12) r = d (cid:19) ∆ + . . . (cid:21) + 4 π (cid:90) Rd r e − βu sf ( r ) dr. (10)In the limit ∆ → → ε → ∞ ,then some of the expansion terms must also be retained.The correct way to carry out the limit is to require that∆ e βε = const, which is referred to as the Baxter limit ,and yieldslim ∆ → ε →∞ Z = 4 πd l g e − βu sf ( d ) + 4 π (cid:90) Rd r e − βu sf ( r ) dr, (11)where we introduced the “sticky length” defined as l g = lim ε →∞ ∆ → ∆ e βε . (12)Of great concern for simulations is the width ∆ of thewell potential, since in practice the exact Baxter limitcannot be attained and ∆ must remain finite. To esti-mate what is sufficiently small value of ∆, we considerthe previous simple system with u sf = 0, for which theexact partition function is Z = 4 πd l g (cid:20) d + 13 (cid:18) ∆ d (cid:19) (cid:21) + 4 πR (cid:20) − (cid:18) d + ∆ R (cid:19) (cid:21) . (13)If we ignore the second term in square brackets, assum-ing R (cid:29) d , we conclude that the well potential becomessticky if ∆ /d (cid:28)
1. In practice, we find that ∆ /d ≈ .
01 issufficiently small to suppress most contributions of finite∆.The Baxter sticky potential may appear analogous to adelta function potential often used in quantum mechan-ics. This, however, is misleading. The well potential inEq. (7) transforms into the delta function in the limits∆ → ε → ∞ , while the product ∆ ε is held fixed.On the other hand, the Baxter limit, requires that ∆ e βε remains constant. To see this more clearly we define f ( r ) = (cid:40) , d ≤ r ≤ d + ∆ , , r < d or r > d + ∆ . (14)The Boltzmann factor then can be written as e − βu well ( r ) = 1 + ∆(e β(cid:15) − f ( r ) , (15)which in the Baxter limit reduces tolim ε →∞ ∆ → e − βu well ( r ) = 1 + l g δ ( r − d ) , (16)with the sticky length given by l g ≡ ∆(e β(cid:15) − − l g can be ne-glected, however, in the simulations with finite ∆ we willuse the exact expression for l g . The sticky potential itselfis then βu well ( r ) = − ln (cid:20) l g δ ( r − d ) (cid:21) , (17)which shows that it is weaker than the delta functionpotential. Indeed, a delta function potential would resultin an irreversible association between the sticky spheres.Finally, we note that if the expression (16) is used in thepartition function Eq. (9), we will arrive directly at theEq. (11). B. Monte Carlo simulations details
We are now in a position to implement numerical simu-lations. The simulations are performed inside a sphericalWigner-Seitz (WS) cell of radius R . A spherical colloidalparticle of radius a is placed in the center of the WScell. The radius of the cell is determined by the colloidalvolume fraction of the suspension, ϕ c = a /R . The mo-tivation for using WS is that for small salt concentrationcolloidal system may crystallize, in which case thermody-namics will be very well described by the WS cell model,with a Donnan potential used to control the charge neu-trality. In fact, even for the disordered state the WSapproach to thermodynamics is found to lead to osmoticpressures in excellent agreement with experiments . Ina sense, the many body colloid-colloid interactions in thegrand-canonical ensemble are all included through theboundary condition of vanishing electric field at the WScell boundary.The colloidal particle has N site adsorption sites ran-domly distributed over its surface. Each adsorption siteis a sphere of diameter d , see Fig. 2. If an active siteis basic – has zero charge — it interacts with the hydro-nium ions through the hard core and the Baxter stickypotential, Eq.(8). On the other had if the site is acidic— has charge − q , where q is the proton charge — in ad-dition to the Baxter and hard core interactions, there isalso a long range Coulomb potential between the adsorp-tion site and all the ions inside simulation cell. In thiswork all the ions and the adsorption sites have diameter4 ˚A and the colloidal particle has radius of a = 100 ˚A.The system is connected to a reservoir of strong acid atconcentration 10 − pH , and a reservoir of 1:1 strong elec-trolyte at concentration c s . The solvent is considered tobe a uniform dielectric of permittivity (cid:15) w = 80 (cid:15) and theBjerrum length is λ B = q /(cid:15) w k B T = 7 . U = (cid:88) i>j q i q j (cid:15) w | r i − r j | + (cid:88) (cid:48) u well ( r i ) , (18)where the first sum is over all the charged particles, in-cluding the adsorption sites, and the second sum is for FIG. 2: The Baxter’s sticky spherical sites on thecolloidal surface.the sticky interaction between the hydronium ions andthe adsorption sites. The hardcore interaction betweenions, sites, and colloidal surface is implicit. The restric-tion on the second sum indicated by the prime is due tothe fact that each functional site can adsorb at most onehydronium ion. This is the case for carboxyl or aminegroups. Therefore, once there is a hydronium ion withinthe distance ∆ of the adsorption site, the short rangesticky potential of this site with other hydronium ions isswitched off. In this paper we will not consider more com-plicated metal oxide ions which can adsorb more than oneproton. To perform simulations we used Metropolis al-gorithm . For large WS cells, when a system establishesa well defined bulk concentration far from the colloidalsurface, we can use canonical Monte Carlo simulations .For large colloidal volume fractions, when WS is smalland bulk concentration is not reached inside the cell, weuse the grand canonical Monte Carlo simulations . Thisis done in order to have a well defined reservoir concen-trations of acid c a and salt c s , which are necessary tocompare the theory with the simulations. In both typesof simulations we have used 5 × MC steps for equili-bration and 10 steps for production.We first check the convergence of MC results to theBaxter sticky limit by studying systems with differentvalues of ∆ and (cid:15) , while keeping fixed the sticky length l g . Fig. 3, shows the rapid convergence to the Baxterlimit, with decreasing value of ∆. III. EQUILIBRIUM CONSTANT FOR PARTICLESINTERACTING VIA A PAIR POTENTIAL
To connect the simulations presented in the previoussection with the NP theory we must relate the stickylength with the bulk association constant.FIG. 3: The density profiles of ions for different size ∆and fixed l g =109.9˚A. The density is plotted in terms ofthe number of particles per ˚A . A. Neutral pairs
We first consider a general two-component system ofsticky spheres. To avoid confusion with previous labels,we designate the “atoms” of each species as X and Y.The interaction between atoms of the same species is u xx ( r ) = u yy ( r ) = u hs ( r ) + u sf ( r ) , (19)and between the atoms of different species is u xy ( r ) = u hs ( r ) + u sf ( r ) + u well ( r ) , (20)This is the, so called, “physical picture”, in which onlyatoms exist. Alternatively, we can regard two atoms X and Y in contact to form a molecule XY . This corre-sponds to the “chemical picture”, see Fig. (4) for illus-tration. In the chemical picture, we have free atoms X and Y , and molecules XY which are in “chemical” equi-librium , X + Y (cid:10) XY . (21)At most two atoms X and Y are permitted to interact viaa sticky potential. Without this restriction, one has to ac-count for the presence of triplets XYX, quartets XYXY,and other higher order formations, together with theircorresponding chemical reactions.To obtain the equilibrium constant for the chemicalreaction in Eq. (21) we compare the equations of statecalculated using the physical and the chemical pictures.Clearly the osmotic pressure calculated using the two in-terpretations of the same physical reality has to be same.Within the physical interpretation, the system is com-prised of two types of atoms, X and Y, and the formationof pairs XY is devoid of any special meaning. The virial FIG. 4: The two representations correspond to a) aphysical and b) a chemical interpretation.expansion of the osmotic pressure up to second order inbulk concentration c i is βP phys = c x + c y + B xx c x + B yy c y +2 B xy c x c y + . . . , (22)where B ij are the second virial coefficients defined as B ij = 2 π (cid:90) ∞ (cid:16) − e − βu ij ( r ) (cid:17) r dr, (23)and whose various contributions are B xx = B yy = B hs + B sf ,B xy = B hs + B st + B sf , (24)which, after evaluation, become B hs = 2 πd ,B st = − πl g d e − βu sf ( d ) ,B sf = 2 π (cid:90) ∞ d (cid:16) − e − βu sf ( r ) (cid:17) r dr,B sf = 2 π (cid:90) ∞ d (cid:16) − e − βu sf ( r ) (cid:17) r dr, (25)where B st was evaluated using the Boltzmann factor inEq. (16). Inserting these contributions into the expan-sion in Eq. (22) yields βP phys = c x + c y + B hs (cid:0) c x + c y (cid:1) + B sf (cid:0) c x + c y (cid:1) + 2 B sf c x c y + 2 B st c x c y + . . . , (26)where the first line is the ideal-gas contribution, the sec-ond line is the second order correction due to hard-coreand soft interactions, and the third line is the secondorder correction due to the sticky potential.To formulate the equation of state within the chemi-cal picture, we need to define the concentrations of freeatoms X, Y, and of molecules XY, designated by the su-perscript *: c ∗ x = c x − c ∗ xy ,c ∗ y = c y − c ∗ xy ,c ∗ xy = K Bulk c ∗ x c ∗ y , (27)where the last equation was obtained using the definitionof the equilibrium constant of the reaction in Eq. (21), K Bulk = c ∗ xy c ∗ x c ∗ y , (28)valid in the dilute limit where activities are approximatedby concentrations. To second order in concentrations c i ,Eq. (27) can be written as c ∗ x = c x − K Bulk c x c y + . . . ,c ∗ y = c y − K Bulk c x c y + . . . ,c ∗ xy = K Bulk c x c y + . . . . (29)In the chemical picture, the interactions between freeatoms X and Y do not include sticky interaction, u (cid:48) xy ( r ) = u hs ( r ) + u sf ( r ) . which acts only within a molecule XY. Without the stickyinteraction, the modified second virial coefficient in thechemical interpretation is B (cid:48) xy = B hs + B sf . (30)The osmotic pressure up to second order in concentra-tions c i is then βP chem = c ∗ x + c ∗ y + c ∗ xy + B xx c ∗ x + B yy c ∗ y +2 B (cid:48) xy c ∗ x c ∗ y + . . . . (31)The terms2 B x,xy c ∗ x c ∗ xy + 2 B y,xy c ∗ y c ∗ xy + B xy,xy c ∗ xy , (32)that are second order in c ∗ i are omitted since, due to c ∗ xy ≈ K Bulk c x c y in Eq. (29), they are of higher order in c i . Using formulas in Eq. (29), and substituting for thecoefficients B ij , Eq. (31) becomes βP chem = c x + c y + B hs (cid:0) c x + c y (cid:1) + B sf (cid:0) c x + c y (cid:1) + 2 B sf c x c y − K Bulk c x c y + . . . , (33)Setting P phys = P chem , and matching the terms of thesame order yields K Bulk = − B st = 4 πl g d e − βu sf ( d ) , (34)We note that the above derivation assumes a dilutelimit, where the definition of K Bulk in Eq. (28) and thesecond order expansion of βP are valid. The result in Eq. (34), however, is exact for any concentration. Thisis because the quantity K Bulk itself is independent of con-centrations. We simply took advantage of this fact andchose the limit where all the expressions are the simplest.If we set u fs = 0 , the equilibrium constant becomes K Bulk = 4 πl g d , (35)which is appropriate for the pair formation between basesand hydronium ions. On the other hand, as we will see inthe following section, Eq. (34) with u sf correspondingto the Coulomb potential will be appropriate for weakacid-hydronium equilibrium constant. B. Charged pairs
In bulk, acid “molecule” dissociates resulting in a hy-dronium ion and a corresponding anion:HA (cid:29) H + + A – , (36)The thermodynamics of bulk electrolytes, even withoutcovalent bonding between the ions, is complicated by thedivergence of the virial expansion due to the long rangenature of the Coulomb interaction. Instead a certainclass of perturbative diagrams must be summed togetherto obtain a finite result . This leads to a non-analyticterm in the density expansion of the osmotic pressurewhich scales with electrolyte concentration as c / . Thenext order term which scales as c can be interpretedas the result of Bjerrum anion-cation pair formation. Inthe case of purely electrostatic interactions, the equilib-rium constant for such cluster formation was derived byEbeling considering the exact density expansion of theequation of state up to O (cid:16) c / (cid:17) . The Ebeling equilib-rium constant is: K Eb = 8 πd (cid:26) b (cid:2) Ei ( b ) − Ei ( − b ) (cid:3) −
13 cosh b − b sinh b − b cosh b + 13 + 12 b (cid:27) , (37)where b = λ B d . For large values of b (strong couplinglimit), the equilibrium constant can be expanded asymp-totically to give: K Eb = 4 πa e b b (cid:18) b + 4 × b + 4 × × b + ... (cid:19) . (38)This may be compared with the Bjerrum phenomeno-logical association constant for formation of anion-cationpairs K Bj = 4 π (cid:90) R Bj d e λBr r dr, (39)where R Bj = λ B / K Bj is completely in-sensitive to the precise value of cutoff R Bj . Further-more, the low temperature expansions for K Eb and K Bj are found to be identical . One can then interpret theEbeling equilibrium constant as the analytic continua-tion of K Bj over the full temperature range. With thisobservation it becomes easy to obtain the equilibriumconstant for sticky electrolytes. In the spirit of Bjerrum,we then write K Bulk = 4 π (cid:90) R Bj d e − βu st ( r )+ λBr r dr. (40)Using Eq.(16) we obtain K Bulk = 4 π (cid:90) R Bj d (cid:2) (1 + l g δ ( r − d ) (cid:3) e λBr r dr , (41)which after integration yields, K Bulk = 4 πd l g e b + (cid:90) R Bj d e λBr r dr . (42)The validity of the above equation is extended beyondthe strong coupling limit by replacing the integral with K Eb . In the case of weak acids, large l g , the first term willdominate Eq. (42), so that the bulk equilibrium constantfor a weak acid can be approximated by K Bulk = 4 πd l g e b , (43)which is similar to Eq. (34) of the previous section. IV. UNIFORMLY STICKY COLLOID
To build a theory of charge regularization of colloidalparticles we start with the simplest possible model inwhich the whole of colloidal surface is sticky. Colloidalparticle of radius a is placed at the center of a sphericalWS cell of radius R . The density profiles of ions, then,satisfy the modified Poisson-Boltzmann (mPB) equation: ∇ φ ( r ) = − π(cid:15) w σ δ ( r − a − r ion ) − πq(cid:15) w (cid:2) c H + ( r ) + c + ( r ) − c − ( r ) (cid:3) , (44)where σ = 0 if the surface groups are basic, and σ = − N sites q/ π ( a + r ion ) if all the groups are acidic. Theionic concentrations are defined as: c H + ( r ) = c a e − β ( u ( r )+ qφ ( r ) ) (45) c + ( r ) = c s e − βqφ ( r ) (46) c − ( r ) = ( c a + c s ) e βqφ ( r ) , (47)where u ( r ) is the sticky potential between the colloidalsurface and a hydronium ion, c a = 10 − pH is the reser-voir concentration of acid, and c s is the reservoir con-centration of 1:1 salt. We assume that both acid andsalt in the reservoir are strong electrolytes and are fullydissociated. Using Eq. 16 we obtain e − β ( u ( r )+ qφ ( r ) ) =(1 − l g δ ( r − a − r ion ))e − βqφ ( r ) , which means that the sur-face density of adsorbed hydronium ions is : σ su = qc a l g e − βqφ , (48) where φ = φ ( a + r ion ). The net surface charge densityis then σ net = σ + σ su . (49)To calculate the ionic density profiles and the numberof condensed hydronium ions we must now solve the PBequation ∇ φ ( r ) = 8 πq(cid:15) w ( c a + c s ) sinh[ βφ ( r )] , (50)with the boundary conditions φ (cid:48) ( R ) = 0 and φ (cid:48) ( a + r ion ) = 4 πσ net /(cid:15) w . The calculation can be performednumerically using the 4th order Runge-Kutta, in whichthe value of the surface potential φ ( a + r ion ) = φ is ad-justed based on the Newton-Raphson algorithm to obtainzero electric field at the cell boundary.In realty, however, the whole of colloidal surface is notuniformly sticky and hydronium ions can only adsorb onspecial sites . We now explicitly consider the modifica-tions that must be made to the above theory in order toaccount for the discrete nature of adsorption sites. V. NEUTRAL FUNCTIONAL GROUPS
We first consider a colloidal particle with N site neu-tral basic groups (sticky spheres) uniformly distributedon its surface. To simplify the geometry we will map thespherical sticky sites onto circular sticky patches of thesame effective contact area. Since both hydronium andthe adsorption sites are modeled by spheres of the samediameter d , the hard core repulsion between the colloidalsurface and the hydronium ion restricts the effective con-tact area to 2 πd . Therefore, the patch radius must be r patch = √ d, (51)Compared to the situation discussed in the previous sec-tion in which the whole of colloidal surface was sticky,the effective area on which hydronium ions can becomeadsorbed is significantly reduced in the case of discreteadsorption sites . Nevertheless, we can still use the sameapproach as in Section IV, if the sticky length is rescaledas l effg = l g α eff , to account for the reduced adsorptionarea, where α eff = N actsite πr patch π ( a + r ion ) , (52)is the fraction of the surface area occupied by the active sticky patches. Note that if hydronium is adsorbed to apatch, this patch becomes inactive, preventing more thanone hydronium ion from being adsorbed. The numberof adsorbed hydronium ion is given by Eq. (48) with l g replaced by l effg . As the process of adsorption progresses,the number of active sites decreases in such a way as N actsite = N site − π ( a + r ion ) c a l effg e − βφ , (53)resulting in a self-consistent equation for l effg . SolvingEqs. (52) and (53), the effective sticky length is found tobe l effg = l g N site r patch a + r ion ) (cid:16) l g c a e − βφ πr patch (cid:17) . (54)The effective surface charge density which must be usedas the boundary condition for PB equation is then σ eff = qc a l effg e − βφ = qK Surf N site c a e − βφ π ( a + r ion ) (cid:0) K Surf c a e − βφ (cid:1) , (55)where K Surf = 2 πl g d = K Bulk /
2, where the bulk asso-ciation constant is the same as in Eq. 35.We stress again that the bulk equilibrium constant K Bulk is exact for the model of sticky hard spheres anddoes not depend on the density of the reactants. Thehigher order terms of the virial expansion, however, willmodify the activity coefficients, so that in the law of massaction the concentrations will have to be replaced by theactivities. Nevertheless, since the PB equation does nottake into account ionic correlations, to remain consistent,the activity coefficients must also be set to unity. Solv-ing Eq. 50, with the boundary conditions φ (cid:48) ( R ) = 0 and φ (cid:48) ( a + r ion ) = 4 πσ net /(cid:15) w we are able to obtain the densityprofile of ions around the colloidal particle.To explore the range of validity of the theory we willcompare it with the results of Monte Carlo simulations.We first consider colloidal particles with 300 and 600 ac-tive neutral basic sites and concentration of HCl set to50 mM. In Figs. 5 and 6 we show the comparison be-tween the simulation data, NP theory, and the presentwork, For these parameters the difference between thenew theory and the NP approach is not very large, nev-ertheless it is clear that the simulation results are in amuch better agreement with the theory developed in thepresent paper. The figures show that the density of freehydronium ions decreases near the colloidal surface. Thisis not surprising, since once some of the hydroniums haveadsorbed to the neutral basic groups, colloidal surfacebecomes positively charged and repels other cations. Amore curious behavior is found for the anion Cl – , theconcentration of which shows a peak close to the surface,but then diminishes on further approach. The reason forthis is that anions prefer to stay close to the adsorbedcations H + , which in turn want to minimize the repul-sive electrostatic energy between themselves, as well asto maximize entropy. This favors the hydronium ions tobe located at about 3 r ion from the colloidal surface.This is precisely the position of the peak found in thedensity profile of anions. This fine detail, however, is be-yond the scope of the present theory. Nevertheless thefact that the density profiles away from colloidal surfaceare perfectly described by the present theory implies thatthe prediction for the total number of adsorbed hydro-nium ions is correct, in spite of the fine structure of ionicdensity profiles near the surface.
100 120 140 160 180 r [Å] × -5 × -5 × -5 × -5 c MC SimulationTheory of current workNP Theory 100 120 140 160 180 r[Å] × -5 × -5 × -5 a) b) FIG. 5: Density profiles of hydronium and Cl – measuredin particles per ˚A . Symbols are the simulation dataand solid (green) and dashed (blue) lines are thepredictions of the NP theory and of the theorydeveloped in the present work, respectively. Theparameters are a = 100 ˚A, R = 200 ˚A, and l g = 109 . – and b) Density profileof hydronium.
100 120 140 160 180 r [Å] × -5 × -5 × -5 × -5 × -5 × -5 × -5 c MC SimulationTheory of current workNP Theory 100 120 140 160 180 r[Å] × -5 × -5 × -5 a) b) FIG. 6: Density profiles of hydronium and Cl – measuredin particles per ˚A . Symbols are the simulation dataand solid (green) and dashed (blue) lines are thepredictions of the NP theory and of the theorydeveloped in the present work, respectively. Theparameters are a = 100 ˚A, R = 200 ˚A, and l g = 109 . – and b) Density profile ofhydronium.We next consider the effect of 1:1 salt on the chargeregulation. We study a colloidal particle with 200 basicsites in the presence of HCl and NaCl, both at concen-tration 10 mM. We assume that both acid and salt arecompletely ionized. The results of the theory and sim-ulations are shown in Fig. 7. Once again we see a verygood agreement between the present theory and the MCs
100 120 140 160 1801 × -5 × -5 c MC simulationTheory o of current workNP Theory100 120 140 160 1805 × -6 × -6 × -6 c
100 120 140 160 180 r [Å] × -6 × -6 × -6 c Cl - H + Na + FIG. 7: Density profiles of hydronium, Cl – , and Na + measured in particles per ˚A . Symbols are thesimulation data and solid (green) and dashed (blue)lines are the predictions of the NP theory and of thetheory developed in the present work, respectively. Theparameters are a = 100 ˚A, R = 200 ˚A, and l g = 109 . .simulations.In experiments, Zeta potential is more easily availablethan the effective charge. Definition of Zeta potential,however, requires knowledge of the position of the slipplain. Nevertheless we expect that Zeta potential willbehave similarly to the electrostatic contact surface po-tential. In Figs. 8 and 9 we show the behavior of thesurface potential and the effective charge Z eff in unit ofcharge q as a function of pH, for the present theory andNP theory, respectively and in Fig 10, the behavior ofthe two as a function of salt concentration. We observethat addition of 1:1 electrolyte diminishes the contactpotential. This, in turn, lowers the electrostatic energypenalty for bringing hydronium ions to colloidal surface,thus favoring their association with the active sites. In-deed, Fig 10b shows that the effective charge of colloidalparticle increases with increasing salt concentration. VI. CHARGED FUNCTIONAL GROUPS
We next consider a colloidal particle with N site acidicsurface groups each carrying a charge − q . If all thegroups would be ionized, the particle would acquire anet charge Q = − N site q . The chemical equilibrium be-tween hydronium and acid groups, however, reduces thisvalue to Q eff = Q + Q con , where Q con = 4 π ( a + r ion ) qc a l effg e − βqϕ (56)is the number of associated hydronium ion and ϕ is thepotential of mean force (PMF) — the work required tobring an ion from the bulk to contact with one of the pH q | φ C on t ac t | / k b T NPPresent Theory N Site = 300
FIG. 8: Contact potential as a function of pH for thepresent theory and NP theory, respectively. Thecolloidal particle has 300 basic functional group on itsurface. The parameters are a = 100 ˚A, R = 200 ˚A, and l g = 109 . pH | Z e ff | NPPresent Theory N site = 300 FIG. 9: Effective charge of colloidal particle in unit of q as a function of pH for the present theory and NPtheory, respectively. The colloidal particle has 300 basicactive functional group on it surface. The parametersare a = 100 ˚A, R = 200 ˚A, and l g = 109 . φ and a contribution fromthe discrete nature of surface charge groups µ qqc , ϕ = φ + µ qqc . (57)The value of l effg is given by Eq. (54) with the mean-fieldpotential replaced by the PMF, φ → ϕ . The effectivesurface charge density then reduces to σ eff = − N site q π ( a + r ion ) + qK aSurf N site c a e − βφ π ( a + r ion ) (cid:16) K aSurf c a e − βφ (cid:17) (58)where K aSurf = K aBulk − b − βµ qqc , (59)0 c s [M] q | φ C on t ac t / k b T | pH = 2pH = 3 c s [M] Z e ff a) b) FIG. 10: (a) Contact potential and (b) the effect chargeof colloidal particle in unit of q as a function of saltconcentration c s (M) for different pH. The colloidalparticle has 300 basic functional group. The parametersare a = 100 ˚A, R = 200 ˚A, and l g = 109 . K aBulk is given byEq. (43). The term e − b in Eq. 59 discounts the directCoulomb interaction between the hydronium ion and itsadsorption site, which is already accounted for in the µ qqc . VII. THE EFFECT OF DISCRETE CHARGES
It is well known that the PB equation is very accuratefor systems containing only 1:1 electrolyte. The mean-field nature of this equation is manifested by the com-plete neglect of ionic correlations, which are found to besmall for aqueous solutions of monovalent ions . How-ever, in the case of acidic groups, hydronium ions willcondense directly onto charged sites and discrete natureof hydronium ions and surface sites can not be neglectedfor the associated ions. The free ions, however, can stillbe treated at the mean-field level.To account for the discrete nature of surface groups,we add and subtract a uniform neutralizing backgroundto the colloidal surface. The negative of the backgroundcan be combined with the mean-field electrostatic poten-tial produced by the ions to yield the total mean-fieldelectrostatic potential φ ( r ). The potential produced bythe discrete surface charge and their neutralizing back-ground, on the other hand, correspond to µ qqc defined inEq.(57). To calculate µ qqc we will ignore the curvature ofthe colloidal surface. Furthermore, we will suppose thatthe adsorption sites are uniformly distributed, forming atriangular lattice of spacing L .We start by calculating the electrostatic potential pro-duced by an infinite planar triangular array of charges,see Fig. 11. This potential must satisfy the Poisson equa-tion ∇ G ( r ) = − πq(cid:15) w (cid:88) n,m δ ( z ) δ ( ρ − n a − m a ) , (60) FIG. 11: The triangular lattice used to evaluate µ qqc and µ qqn .where z and ρ = x ˆ x + y ˆ y are the transverse and lon-gitudinal directions, respectively, and the lattice vectorsare given by a = L ˆ x , a = 12 L ˆ x + √ L ˆ y . (61)The area of the unit cell of triangular lattice is | γ | = | a × a | = √ L (62)The reciprocal lattice vectors b i are defined as a j · b j =2 πδ ij , and are given by b = 2 πd (cid:18) ˆ x − ˆ y √ (cid:19) . b = 2 πd (cid:18) ˆ y √ (cid:19) . (63)The periodic delta function can be written as (cid:88) n,m δ ( ρ − n a − m a ) =1 γ (cid:88) n,m e i b · ρ n + i b · ρ m (64)and the Green function as G ( r ) = 1 γ (cid:88) n,m g n,m ( z )e i b · ρ n + i b · ρ m , (65)where g n,m ( z ) is a function of z coordinate only. Substi-1tuting Eq. (65) into Eq. (60) we obtain ∂ g n,m ( z ) ∂z − k g n,m ( z ) = − πq(cid:15) w δ ( z ) ,k = (cid:118)(cid:117)(cid:117)(cid:116) π L (cid:32) n + (cid:18) m √ − n √ (cid:19) (cid:33) , (66)which has a solution of the form g n,m ( z ) = (cid:40) A e − kz , z > ,A e kz , z < , (67)Integrating Eq. (66) once, we see that the derivative of g ( z ) is discontinuous at z = 0 with g (cid:48) n,m (0 + ) − g (cid:48) n,m (0 − ) = − πq(cid:15) w , (68)from which we determine A = 2 πq/(cid:15) w k ,. The Greenfunction can then be written as G ( r ) = 2 πqγ(cid:15) w n = ∞ (cid:88) n = −∞ m = ∞ (cid:88) m = −∞ e − k | z | k cos 2 πL (cid:18) nx + 1 √ ym − yn ) (cid:19) . (69)The ( n = 0 , m = 0) term of G ( r ) diverges. Indeed, if wetake the limit k → | z | . This is nothing more than the potential ofa uniformly charged plane. Therefore, if we introducea neutralizing background, we will cancel precisely thisterm, eliminating the divergence. The electrostatic po-tential produced by a triangular array of charges on aneutralizing background is then¯ G ( r ) = 2 πqγ(cid:15) w n = ∞(cid:48) (cid:88) n = −∞ m = ∞(cid:48) (cid:88) m = −∞ e − k | z | k cos 2 πL (cid:18) nx + 1 √ ym − yn ) (cid:19) , (70)where the prime on the sums indicates that we have re-moved the term ( n = 0 , m = 0). Bringing an ion ofopposite charge into contact with one of the adsorptionsites then yield µ qqc = − πq γ(cid:15) w n = ∞(cid:48) (cid:88) n = −∞ m = ∞(cid:48) (cid:88) m = −∞ e − kr ion k . (71)Even if sites are not perfectly ordered on the colloidalsurface, we still expect that µ qqc derived in Eq. (71)will provide a reasonably accurate account of the dis-creteness effects assuming that the average separationbetween Z acid sites is such that the area per site is γ = 4 π ( a + r ion ) /Z , where γ is given by Eq. 62. The
100 120 140 160 180 200 r [Å] -5.0 -4.7 -4.5 -4.4 -4.3 -4.2 -4.2 -4.1 -4.0 c NP TheoryPresented TheoryMC Simulation100 120 140 160 180 200 r [Å] -5.7 -5.4 -5.2 -5.1 c Z = 600 pH = 2 H + Cl - FIG. 12: Comparison between the present theory (solidlines), NP theory (dashed lines) and simulations(symbols), for colloidal particles with Z = 600 and l g = 109 .
100 150 200 250 300 -5.3 -5.0 -4.8 -4.7 c NP TheoryPresented thoeryMC simulation
100 150 200 250 3000.010 -5.0 -4.7 -4.5 -4.4 -4.3 c
100 150 200 250 300 r [Å] -5.7 -5.4 -5.2 -5.1 c Z= 300 pH = 3 H + Na + Cl - (a)(b)(c) FIG. 13: Comparison between the present theory (solidlines), NP theory (dashed lines) and simulations(symbols), for colloidal particles with Z = 300 chargedfunctional group with l g = 109 . .average separation between acid groups is then L =( a + r ion ) (cid:113) π/ √ Z We first consider a colloidal particles with 600 acid sur-face groups with l g = 109 .
97 ˚A, in a solution of pH = 2.The ionic density profiles are presented in Fig. 12. Wesee that the theory is in excellent agreement with simu-lations, while NP approach shows significant deviation.Next we consider particles with 300 charged sites insidean acid solution containing 1:1 salt. Once gain there isa good agreement between theory and simulations, seeFig. 13.In Fig. 14 we show the behavior of the effective chargeand contact potential of colloidal particle as a function of1:1 salt concentration for different pH values. The figure2 q | φ C on t ac t | / k B T pH = 3pH = 2pH = 10 0.05 0.1 C s [M] | Z e ff | FIG. 14: Modulus of the effective charge in unit of q and contact potential of colloidal particle as a functionof 1:1 salt concentration C s for different values of pH.The colloidal particle has 300 charged functional groupon it surface. The parameters are a = 100 ˚A, R = 200˚A, and l g = 109 . shows that increase of salt concentration leads to increaseof the modulus of the effective charge. This, again, is aconsequence of electrostatic screening produced by salton the Coulomb interaction between hydronium ions andthe negatively charged adsorption sites — making theassociation of a hydronium with an active site less ener-getically favorable. In Fig. 15 we compare the effectivecharge and contact potential calculated using the presenttheory and the values predicted by the NP theory, fornanoparticles with 300 charged groups. As can be seen,neglect of discrete charge effects in the NP theory leadsto smaller modulus of the contact potential and of theeffective charge. We also note that at large pH the effec-tive charge saturates at the value smaller than the barecharge. This is a consequence of the overall charge neu-trality of the colloidal suspension. Even if the reservoirhas a very small concentration of acid – large pH, inthe absence of other cations inside the suspension, theremust be enough hydronium ions to compensate all thecolloidal charge. Some of these hydronium ions will thenassociate with the surface groups, leading to the satura-tion of the effective colloidal charge. We now performthe same calculation, but in the present of a reservoirwith 10 mM monovalent salt. As can be seen in Fig. 16,in the presence of salt, for high pH both NP and ourtheory predict that the effective charge approaches thebare charge. This is should be contrasted with the no-salt system. When the system is connected to both thesalt and acid reservoirs, at large pH the hydronium ionsinside the system are replaced by the salt cations, whichthen control the overall charge neutrality of the colloidalsuspension. Since in our model salt cations do not reactwith the surface groups, for reservoir at large pH veryfew hydronium ions will be present inside the suspension.Therefore, all the surface groups will become ionized, and pH q | φ c on t ac t | / k b T The Present TheoryNP 1 1.5 2 2.5 3 3.5 4 pH | Z e ff | FIG. 15: The modulus of the effective charge in units of q and the contact potential of a nanoparticle as afunction of pH in the acid reservoir, predicted by theNP and the present theories. The colloidal particle has300 charged functional group on it surface. Theparameters are a = 100 ˚A, R = 200 ˚A, and l g = 109 . pH q | φ c on t ac t | / k b T NPThe Present Theory 3 3.5 4 4.5 5 5.5 pH | Z e ff | FIG. 16: The modulus of the effective charge in unitsof q and the contact potential of a nanoparticle as afunction of pH in the acid reservoir, predicted by theNP and the present theories. The colloidal particle has300 charged functional group on it surface. Thesuspension is in a contact with a monovalent saltreservoir at concentration of 10 mM. The parametersare a = 100 ˚A, R = 200 ˚A, and l g = 109 . VIII. MIXTURE OF FUNCTIONAL GROUPS
As a final example we consider a colloidal particle witha mixture of acidic and basic surface groups. Followingthe same approach introduced in the previous sections3we find that the effective surface charge is σ eff = − N acid q π ( a + r ion ) + q ( l effgc c a e − βϕ c + l effgn c a e − βϕ n ) ,l effgc = l gc N acid r patch a + r ion ) (cid:16) l gc c a e − βϕ c πr patch (cid:17) ,l effgn = l gn N base r patch a + r ion ) (cid:16) l gn c a e − βϕ n πr patch (cid:17) , (72)where N acid is the number of acidic groups and N base = N site − N acid is the number of basic groups. The effec-tive sticky length for acidic (charged) and basic (neutral)groups are: l effgc and l effgn , respectively. The discretenesseffects will manifest themselves in different ways for hy-dronium ions condensing on acidic and basic groups, βϕ c = βφ c + µ qqc ,βϕ n = βφ n + µ qqn , (73)Since the values µ qqc,n depend only on the electrostaticinteraction between the hydronium ion and the charged(acid) sites, the value of µ qqc will be the same as in Eq.(71), depending only on the average separation betweenthe acidic groups. We will suppose that the basic groupsare also uniformly distributed on the colloidal surfaceon a dual hexagonal lattice with vertexes at the centerof each triangle composed of acidic sites. In this casethe position of one of the basic groups will be at x = d/ y = √ d/
4. Using Eq. (70) we obtain µ qqn µ qqn = − πq γ(cid:15) w n = ∞(cid:48) (cid:88) n = −∞ m = ∞(cid:48) (cid:88) m = −∞ e − kr ion k cos π (cid:18) m + n (cid:19) (74)The effective surface charge density can now be writtenas σ eff = − N acid q π ( a + r ion ) + qK aSurf N acid c a e − βφ π ( a + r ion ) (cid:16) K aSurf c a e − βφ (cid:17) + qK bSurf N base c a e − βφ π ( a + r ion ) (cid:16) K bSurf c a e − βφ (cid:17) (75)where K aSurf = K aBulk − b − βµ qqc , (76)and K bSurf = K bBulk − βµ qqn . (77)The bulk equilibrium constants for acid and base, K aBulk and K bBulk , are given by Eqs. (43) and (35), respectively.To test our theory for mixture of basic and acidic surface N acid
500 450 300 200 50 µ qqc -0.6278 -0.6653 -0.8073 -0.94227 -1.31329 µ qqn ——— 0.14311 0.1519 0.1529 0.119273 TABLE I: Different values of µ qqc and µ qqn for mixture ofcharged sites. The total number of adsorption sites is500groups we, once again, compare it with MC simulations.We consider a colloidal particle with 500 adsorption siteswith different number of acidic groups N acid . The stickylengths — l gc and l gn — are 109 .
97 and 1099 .
7, respec-tively. These values correspond to the equilibrium con-stants K eq = 0 .
012 and 0 . l g increases, it becomes progressively more difficult toequilibrate the simulations. For this reason we have cho-sen values of l g that are not too large. This, however,has no implication for the theory, which remains validfor arbitrary values of l g and K eq .Since the theory is completely general, the values of l g are arbitrarily and one can, in practice, choose thedepth and the width of the sticky potential and calcu-lated the sticky length. There is, however, an additionalconstraint. The equilibration of the simulations becomesprogressively more difficult with increase of sticky length.To have a good test of the theory we, therefore, need tochose sufficiently large sticky length to have a significantassociation of hydronium ions with the adsorption sites,while keeping a reasonable equilibration CPU time. Fur-thermore, to better test the validity of the theory, weshould choose very different sticky lengths for acid andbase sites. This is the reason for a factor of 10 differencebetween the values of l g of acidic and basic groups. In Ta-ble. I we show the values of µ qqc and µ qqn , calculated usingEqs. (71) and (74), respectively – for a colloidal particlewith the total of N site = 500 adsorption sites, N acid ofwhich are acidic (charged) and the rest are basic.Fig. 17 shows that for small number of basic sites, thetheory remains very accurate. This is also the case ifthe number of basic sites is significantly larger than thenumber of acidic sites. The worst agreement is foundwhen N acid ≈ N base in which case the surface of colloidbecomes strongly heterogeneous, with positive, negative,and neutral domains present, leading to the breakdownof assumptions used to calculate µ qqc,n .From the obtained results, we conclude that the the-ory works very well if colloidal particle has either basicor acidic adsorption sites. The discrete charge effects areembedded in the µ qqc and µ qqn , which are calculated usinga regular arrangement of adsorption sites, even thoughin the simulations the sites are randomly distributed.If the number of acid and base sites is approximatelyequal, then after the adsorption, we will end up withlarge domains composed of − µ qqc and µ qqn will break down.4 -4.7 -4.4 -4.2 -4.1 c Cl - H +
100 120 140 160 180 2000.010 -5.0 -4.7 -4.5 -4.4 -4.3
100 120 140 160 180 2000.010 -5.3 -5.0 -4.8 -4.7
100 120 140 160 180 200 r [Å] -5.4 -5.3 -5.2 -5.2 -5.1 -5.0 c
100 120 140 160 180 200 r [Å] -5.3 -5.0 -4.8 -4.7
100 120 140 160 180 200 r [Å] N acid = 500 N acid = 450 N acid = 300N acid = 200 N acid = 50 N acid = 0 FIG. 17: Density profile of ions around colloidal particlewith different number of charged and neutral functionalgroups. The total number of sites is N site = 500, the ofthe reservoir is pH = 2 and there is no additional salt.Symbols are the results of MC simulations and the linesare predictions of the theory. The densities are in unitsof particles per ˚A Nevertheless the theory is found to work quite well, aslong as the number of acidic and basic sites is not thesame. Addition of salt to the system results in even bet-ter agreement between theory and simulations. There-fore, the salt free case, provides the most stringent testof the theoretical approach.
IX. CONCLUSION
In this work, we used a sticky sphere model to mimicchemical reaction on a colloidal surface. Within our the-ory, the discrete charge effects come only from acid sur-face sites, since ions are treated at the mean field Poisson-Boltzmann level. This is the reason why the surface equi-librium constant is dependent only on the value of µ qq .In the previous work we studied colloidal particles withonly acidic surface groups . In that approach we usedone component plasma model (OCP), to account for theelectrostatic corrections due to discrete surface groups.While the approach in Ref. was sufficiently accuratefor colloidal particles with only acidic groups, the OCPdoes not take into account ionic radius, which preventedus from extending this approach to colloidal particleswith a mixture of acidic and basic surface groups. Inthe present study, we have developed a completely differ-ent method to account for discrete surface effects usingperiodic Green functions. The fact that theory workswell even for very small nanoparticles of 10 nm radiusshows the robustness of our approach.The microscopic model presented in the paper permitsus to study the same chemical reaction taking place inbulk and at the interface. In the case of neutral basic surface groups our theory reduces to the NP approachwith the bulk equilibrium constant replaced by the sur-face equilibrium constant, K bSurf → K bBulk . (78)The difference between surface and bulk equilibrium con-stants is a consequence of steric repulsion, which restrictsthe overall surface area of the adsorption sites availablefor interaction with hydronium ions.For colloidal particles with acidic surface groups thesituation is significantly more complex. In this case ourtheory reduces to the NP approach with an effective sur-face equilibrium constant only for weak acidic groups.For such systems we find the surface equilibrium con-stant to be K aSurf = K aBulk − b − βµ qqc , (79)where µ qqc accounts for the discreteness of surface charge.For colloidal particles with a mixture of acidic and basicsurface groups, the respective surface equilibrium con-stants are given by Eqs. (76) and (77).It is important to stress that Eqs. (78) and (79) are notuniversal, and in general will depend on the details of thesystem. These details may include molecular geometry,modified electronic structure of surface functional groups,water structure, etc. Nevertheless the model of stickyadsorption sites demonstrates that there is a mappingbetween the bulk and the surface equilibrium constantswhich allows one to use the Poisson-Boltzmann frame-work to accurately account for the charge regulation incolloidal systems. Any deviations from experiment cantherefore be attributed to the shortfall of the model andnot to the theoretical method used to solve it.In this work our primary goal was to explore the extentof validity of the mean-field NP approach by applyingit to an exactly solvable model. Clearly the microscopicmodel that we used for spherically symmetric hydronium,uniform dielectric water, sticky interactions for covalentbinding, etc., is a very rough approximation to the physi-cal reality. The advantage is that we can solve this modelexactly using computer simulations. Applying the NPapproach to the same model we can then test the extentof validity of the mean-field approximations. We shouldstress that the NP theory does not give us any informa-tion whatsoever about the surface equilibrium constantand assumes it to be the same as the bulk associationconstant. We find, on the there hand, that sticky in-teractions result in a breakdown of the mean-field ap-proximations. Surprisingly, however, we find that all thediscreteness effects can be included in a renormalized sur-face association constant, which our theory predicts ex-plicitly. For our microscopic model the correlations andsteric effects lead to lower surface association constant K Surf , compared to the bulk association constant forthe same acid or base, K Bulk . This means that fewer hy-dronium ions will bind to surface groups, implying that5surface p K a will be smaller than bulk p K a . Within thepresent model, there are two contributions which accountfor the decrease of the association constant at the surface.First, is the steric repulsion from the colloidal surface,which diminishes the access of hydronium to acid andbase groups. Within our model the accessible area forthe charge transfer reaction is lowered by a factor of two,which accounts for the factor of 1 / K Surf , that one needs to use inthe NP theory is actually larger than K Bulk . This meansthat the surface p K a is larger than the p K a of the bulkacid . Since our model already takes into account allthe steric and electrostatic effects at the dielectric con-tinuum level, we must conclude that in order to accountfor the experimental results we must included additionaleffects into the model, such dielectric discontinuity acrossthe colloidal structure, water ordering, quantum natureof proton transfer, etc. The approach that we have devel-oped should allow us to explore these additional effectsin order to understand the mechanisms that lead to theincrease of p K a at colloidal surface. This will be thesubject of the future work. X. ACKNOWLEDGMENT
This work was partially supported by the CNPq. C. Rønne, L. Thrane, P.-O. ˚Astrand, A. Wallqvist, K. V.Mikkelsen, and S. R. Keiding, “Investigation of the temperaturedependence of dielectric relaxation in liquid water by thz re-flection spectroscopy and molecular dynamics simulation,”
TheJournal of Chemical Physics , vol. 107, no. 14, pp. 5319–5331,1997. D. Andelman, “Chapter 12 - electrostatic properties of mem-branes: The poisson-boltzmann theory,” in
Structure and Dy-namics of Membranes (R. Lipowsky and E. Sackmann, eds.),vol. 1 of
Handbook of Biological Physics , pp. 603 – 642, North-Holland, 1995. R. Messina, “Electrostatics in soft matter,”
Journal of Physics:Condensed Matter , vol. 21, p. 113102, feb 2009. A. Abrashkin, D. Andelman, and H. Orland, “Dipolar poisson-boltzmann equation: Ions and dipoles close to charge inter-faces,”
Phys. Rev. Lett. , vol. 99, p. 077801, Aug 2007. Y. Levin, “Electrostatic correlations: from plasma to biology,”
Reports on progress in physics , vol. 65, no. 11, p. 1577, 2002. K. Shen and Z.-G. Wang, “Electrostatic correlations and thepolyelectrolyte self energy,”
The Journal of Chemical Physics ,vol. 146, no. 8, p. 084901, 2017. A. M. Smith, A. A. Lee, and S. Perkin, “The electrostatic screen-ing length in concentrated electrolytes increases with concentra-tion,”
The journal of physical chemistry letters , vol. 7, no. 12,pp. 2157–2163, 2016. R. M. Adar, T. Markovich, A. Levy, H. Orland, and D. Andel-man, “Dielectric constant of ionic solutions: Combined effectsof correlations and excluded volume,”
The Journal of ChemicalPhysics , vol. 149, no. 5, p. 054504, 2018. M. I. Bespalova, S. Mahanta, and M. Krishnan, “Single-moleculetrapping and measurement in solution,”
Current Opinion in Chemical Biology , vol. 51, pp. 113 – 121, 2019. Chemical Ge-netics and Epigenetics Molecular Imaging. D. Frydel, S. Dietrich, and M. Oettel, “Charge renormalizationfor effective interactions of colloids at water interfaces,”
Phys.Rev. Lett. , vol. 99, p. 118302, 2007. D. A. Walker, B. Kowalczyk, M. O. de La Cruz, and B. A. Grzy-bowski, “Electrostatics at the nanoscale,”
Nanoscale , vol. 3,no. 4, pp. 1316–1344, 2011. S. Alexander, P. Chaikin, P. Grant, G. Morales, P. Pincus, andD. Hone, “Charge renormalization, osmotic pressure, and bulkmodulus of colloidal crystals: Theory,”
The Journal of ChemicalPhysics , vol. 80, no. 11, pp. 5776–5781, 1984. E. Trizac, L. Bocquet, and M. Aubouy, “Simple approach forcharge renormalization in highly charged macroions,”
Physicalreview letters , vol. 89, no. 24, p. 248301, 2002. A. Adamson and A. Gast,
Chemistry of Surfaces . John Wiley& Sons, New York, NY, USA,, 1997. D. Wang, R. J. Nap, I. Lagzi, B. Kowalczyk, S. Han, B. A. Grzy-bowski, and I. Szleifer, “How and why nanoparticles curvatureregulates the apparent pka of the coating ligands,”
Journal ofthe American Chemical Society , vol. 133, no. 7, pp. 2192–2197,2011. A. Bakhshandeh, D. Frydel, A. Diehl, and Y. Levin, “Chargeregulation of colloidal particles: Theory and simulations,”
Phys-ical Review Letters , vol. 123, no. 20, p. 208004, 2019. R. Podgornik, “General theory of charge regulation and sur-face differential capacitance,”
The Journal of Chemical Physics ,vol. 149, no. 10, p. 104701, 2018. D. Frydel, “General theory of charge regulation within thepoisson-boltzmann framework: Study of a sticky-charged wallmodel,”
The Journal of Chemical Physics , vol. 150, no. 19,p. 194901, 2019. Y. Avni, D. Andelman, and R. Podgornik, “Charge regulationwith fixed and mobile charged macromolecules,”
Curr OpinElectrochem , vol. 13, pp. 70–77, 2019. K. Linderstrøm-Lang, “On the ionization of proteins,”
CR Trav.Lab. Carlsberg , vol. 15, no. 7, pp. 1–29, 1924. A. Majee, M. Bier, R. Blossey, and R. Podgornik, “Charge reg-ulation radically modifies electrostatics in membrane stacks,”
Physical Review E , vol. 100, no. 5, p. 050601, 2019. Y. Avni, T. Markovich, R. Podgornik, and D. Andelman,“Charge regulating macro-ions in salt solutions: screening prop-erties and electrostatic interactions,”
Soft Matter , vol. 14,pp. 6058–6069, 2018. D. Chan, J. PerTam, L. White, and T. Healy, “Regulation ofsurface potential at amphoteric surfaces during partiele-particleinteraction,”
J. Chem. Soc. Faraday Trans , vol. 71, pp. 1046–1057, 1975. R. Pericet-Camara, G. Papastavrou, S. H. Behrens, andM. Borkovec, “Interaction between charged surfaces on thepoisson- boltzmann level: The constant regulation approxima-tion,”
The Journal of Physical Chemistry B , vol. 108, no. 50,pp. 19467–19475, 2004. H. von Grnberg, “Chemical charge regulation and charge renor-malization in concentrated colloidal suspensions,”
Journal ofColloid and Interface Science , vol. 219, no. 2, pp. 339 – 344,1999. H. G. Ozcelik and M. Barisik, “Electric charge of nanopatternedsilica surfaces,”
Physical Chemistry Chemical Physics , vol. 21,no. 14, pp. 7576–7587, 2019. A. Loˇsdorfer Boˇziˇc and R. Podgornik, “Anomalous multipole ex-pansion: Charge regulation of patchy inhomogeneously chargedspherical particles,”
The Journal of Chemical Physics , vol. 149,no. 16, p. 163307, 2018. F. B. van Swol and D. N. Petsev, “Solution structure effects onthe properties of electric double layers with surface charge regu-lation assessed by density functional theory,”
Langmuir , vol. 34,no. 46, pp. 13808–13820, 2018. T. Markovich, D. Andelman, and R. Podgornik, “Charge regu-lation: A generalized boundary condition?,”
EPL (Europhysics Letters) , vol. 113, p. 26004, jan 2016. M. Krishnan, “A simple model for electrical charge in globu-lar macromolecules and linear polyelectrolytes in solution,”
TheJournal of Chemical Physics , vol. 146, no. 20, p. 205101, 2017. M. Krishnan, “Electrostatic free energy for a confined nanoscaleobject in a fluid,”
The Journal of Chemical Physics , vol. 138,no. 11, p. 114906, 2013. R. A. Hartvig, M. van de Weert, J. Østergaard, L. Jorgensen,and H. Jensen, “Protein adsorption at charged surfaces: Therole of electrostatic interactions and interfacial charge regula-tion,”
Langmuir , vol. 27, no. 6, pp. 2634–2643, 2011. D. C. Prieve and E. Ruckenstein, “The surface potential ofand double-layer interaction force between surfaces character-ized by multiple ionizable groups,”
Journal of Theoretical biol-ogy , vol. 56, no. 1, pp. 205–228, 1976. S. L. Carnie and D. Y. Chan, “Interaction free energy betweenplates with charge regulation: a linearized model,”
Journal ofcolloid and interface science , vol. 161, no. 1, pp. 260–264, 1993. R. Netz, “Charge regulation of weak polyelectrolytes at low-andhigh-dielectric-constant substrates,”
Journal of Physics: Con-densed Matter , vol. 15, no. 1, p. S239, 2002. A. Majee, M. Bier, and R. Podgornik, “Spontaneous symme-try breaking of charge-regulated surfaces,”
Soft Matter , vol. 14,pp. 985–991, 2018. J. E. Hallett, D. A. Gillespie, R. M. Richardson, and P. Bartlett,“Charge regulation of nonpolar colloids,”
Soft matter , vol. 14,no. 3, pp. 331–343, 2018. M. Heinen, T. Palberg, and H. L¨owen, “Coupling between bulk-and surface chemistry in suspensions of charged colloids,”
TheJournal of Chemical Physics , vol. 140, no. 12, p. 124904, 2014. C.-Y. Leung, L. C. Palmer, S. Kewalramani, B. Qiao, S. I.Stupp, M. Olvera de la Cruz, and M. J. Bedzyk, “Crys-talline polymorphism induced by charge regulation in ionicmembranes,”
Proceedings of the National Academy of Sciences ,vol. 110, no. 41, pp. 16309–16314, 2013. G. S. Longo, M. Olvera de la Cruz, and I. Szleifer, “Moleculartheory of weak polyelectrolyte gels: The role of ph and saltconcentration,”
Macromolecules , vol. 44, no. 1, pp. 147–158,2011. R. Zandi, B. Dragnea, A. Travesset, and R. Podgornik, “Onvirus growth and form,”
Physics Reports , 2020. D. Roshal, O. Konevtsova, A. L. Boˇziˇc, R. Podgornik, andS. Rochal, “ph-induced morphological changes of proteinaceousviral shells,”
Scientific reports , vol. 9, no. 1, p. 5341, 2019. L. Javidpour, A. L. Boˇziˇc, R. Podgornik, and A. Naji, “Role ofmetallic core for the stability of virus-like particles in stronglycoupled electrostatics,”
Scientific reports , vol. 9, no. 1, p. 3884,2019. C. T. Zahler and B. F. Shaw, “What are we missing by notmeasuring the net charge of proteins?,”
Chemistry A EuropeanJournal , vol. 25, no. 32, pp. 7581–7590, 2019. A. L. Boˇziˇc and R. Podgornik, “ph dependence of charge multi-pole moments in proteins,”
Biophysical Journal , vol. 113, no. 7,pp. 1454 – 1465, 2017. M. Lund, T. ˚Akesson, and B. J¨onsson, “Enhanced protein ad-sorption due to charge regulation,”
Langmuir , vol. 21, no. 18,pp. 8385–8388, 2005. H. Shen and D. D. Frey, “Effect of charge regulation on stericmass-action equilibrium for the ion-exchange adsorption of pro-teins,”
Journal of Chromatography A , vol. 1079, no. 1, pp. 92 –104, 2005. P. M. Biesheuvel and A. Wittemann, “A modified box modelincluding charge regulation for protein adsorption in a sphericalpolyelectrolyte brush,”
The Journal of Physical Chemistry B ,vol. 109, no. 9, pp. 4209–4214, 2005. N. S. Pujar and A. L. Zydney, “Charge regulation and electro-static interactions for a spherical particle in a cylindrical pore,”
Journal of Colloid and Interface Science , vol. 192, no. 2, pp. 338– 349, 1997. Y. Burak and R. R. Netz, “Charge regulation of interactingweak polyelectrolytes,”
The Journal of Physical Chemistry B ,vol. 108, no. 15, pp. 4840–4849, 2004. R. Kumar, B. G. Sumpter, and S. M. Kilbey, “Charge regulationand local dielectric function in planar polyelectrolyte brushes,”
The Journal of Chemical Physics , vol. 136, no. 23, p. 234901,2012. C. Fleck, R. R. Netz, and H. H. von Grnberg, “Poissonboltz-mann theory for membranes with mobile charged lipids and theph-dependent interaction of a dna molecule with a membrane,”
Biophysical Journal , vol. 82, no. 1, pp. 76 – 92, 2002. M. Lund and B. J¨onsson, “Charge regulation in biomolecularsolution,”
Quarterly Reviews of Biophysics , vol. 46, pp. 265–281, aug 2013. M. L. Grant, “Nonuniform charge effects in protein-proteininteractions,”
Journal of Physical Chemistry B , vol. 105,pp. 2858–2863, apr 2001. A. H. Elcock and J. A. McCammon, “Calculation of weakprotein-protein interactions: The pH dependence of the secondvirial coefficient,”
Biophysical Journal , vol. 80, pp. 613–625, feb2001. M. Lund and B. J¨onsson, “On the Charge Regulation of Pro-teins,”
Biochemistry , vol. 44, no. 15, pp. 5722–5727, 2005. A. C. Mason and J. H. Jensen, “Protein-protein binding is oftenassociated with changes in protonation state,”
Proteins: Struc-ture, Function, and Bioinformatics , vol. 71, pp. 81–91, apr 2008. B. Aguilar, R. Anandakrishnan, J. Z. Ruscio, and A. V.Onufriev, “Statistics and physical origins of pK and ionizationstate changes upon protein-ligand binding,”
Biophysical Jour-nal , vol. 98, pp. 872–880, mar 2010. J. St˚ahlberg and B. J¨onsson, “Influence off charge regulation inelectrostatic interaction chromatography of proteins,”
Analyti-cal Chemistry , vol. 68, no. 9, pp. 1536–1544, 1996. H. K. Tsao, “Electrostatic interaction of an assemblage ofcharges with a charged surface: the charge-regulation effect,”
Langmuir , vol. 16, pp. 7200–7209, sep 2000. A. Kubincov´a, P. H. H¨unenberger, and M. Krishnan, “Interfacialsolvation can explain attraction between like-charged objects inaqueous solution,”
The Journal of Chemical Physics , vol. 152,no. 10, p. 104713, 2020. T. Markovich, D. Andelman, and R. Podgornik, “Complex fluidswith mobile charge-regulating macro-ions,”
EPL (EurophysicsLetters) , vol. 120, p. 26001, oct 2017. R. Podgornik and V. Parsegian, “An electrostatic-surface sta-bility interpretation of the hydrophobic force inferred to occurbetween mica plates in solutions of soluble surfactants,”
Chem-ical physics , vol. 154, no. 3, pp. 477–483, 1991. D. Leckband and J. Israelachvili, “Intermolecular forces in bi-ology,”
Quarterly reviews of biophysics , vol. 34, no. 2, pp. 105–267, 2001. G. Trefalt, S. H. Behrens, and M. Borkovec, “Charge regula-tion in the electrical double layer: Ion adsorption and surfaceinteractions,”
Langmuir , vol. 32, no. 2, pp. 380–400, 2016. P. M. Biesheuvel, M. van der Veen, and W. Norde, “A mod-ified poisson-boltzmann model including charge regulation forthe adsorption of ionizable polyelectrolytes to charged inter-faces, applied to lysozyme adsorption on silica,”
The Journalof Physical Chemistry B , vol. 109, no. 9, pp. 4172–4180, 2005. J.-P. Hsu and B.-T. Liu, “Stability of colloidal dispersions:Charge regulation/adsorption model,”
Langmuir , vol. 15,no. 16, pp. 5219–5226, 1999. J. K. Wolterink, L. Koopal, M. C. Stuart, and W. V. Riems-dijk, “Surface charge regulation upon polyelectrolyte adsorp-tion, hematite, polystyrene sulfonate, surface charge regulation:Theoretical calculations and hematite-poly(styrene sulfonate)system,”
Colloids and Surfaces A: Physicochemical and Engi-neering Aspects , vol. 291, no. 1, pp. 13 – 23, 2006. I. Popa, G. Papastavrou, and M. Borkovec, “Charge regu-lation effects on electrostatic patch-charge attraction inducedby adsorbed dendrimers,”
Phys. Chem. Chem. Phys. , vol. 12, pp. 4863–4871, 2010. P. Gong, J. Genzer, and I. Szleifer, “Phase behavior and chargeregulation of weak polyelectrolyte grafted layers,”
Phys. Rev.Lett. , vol. 98, p. 018302, Jan 2007. S. H. Behrens and M. Borkovec, “Exact Poisson-Boltzmann so-lution for the interaction of dissimilar charge-regulating sur-faces,”
Phys. Rev. E , vol. 60, pp. 7040–7048, Dec 1999. J. G. Kirkwood and J. B. Shumaker, “The influence of dipolemoment fluctuations on the dielectric increment of proteins insolution,”
Proceedings of the National Academy of Sciences ofthe United States of America , vol. 38, no. 10, p. 855, 1952. R. Marcus, “Calculation of thermodynamic properties of poly-electrolytes,”
The Journal of Chemical Physics , vol. 23, no. 6,pp. 1057–1068, 1955. S. Lifson, “Potentiometric titration, association phenomena,and interaction of neighboring groups in polyelectrolytes,”
TheJournal of Chemical Physics , vol. 26, no. 4, pp. 727–734, 1957. N. Adˇzi´c and R. Podgornik, “Field-theoretic description ofcharge regulation interaction,”
The European Physical JournalE , vol. 37, p. 49, Jun 2014. A. M. Smith, M. Borkovec, and G. Trefalt, “Forces between solidsurfaces in aqueous electrolyte solutions,”
Advances in Colloidand Interface Science , vol. 275, p. 102078, 2020. C. Safinya and J. Radler,
Handbook of Lipid Membranes:Molecular, Functional, and Materials Aspects . Taylor & Fran-cis, 2014. B. W. Ninham and V. A. Parsegian, “Electrostatic potentialbetween surfaces bearing ionizable groups in ionic equilibriumwith physiologic saline solution,”
Journal of Theoretical biology ,vol. 31, no. 3, pp. 405–428, 1971. D. L. Chapman, “Li. a contribution to the theory of electro-capillarity,”
The London, Edinburgh, and Dublin philosophicalmagazine and journal of science , vol. 25, no. 148, pp. 475–481,1913. M. Gouy, “Sur la constitution de la charge ´electrique `a la surfaced’un ´electrolyte,”
Anniue Physique(Paris) , vol. 9, pp. 457–468,1910. R. Baxter, “Percus-Yevick equation for hard spheres with sur-face adhesion,”
The Journal of Chemical Physics , vol. 49, no. 6,pp. 2770–2774, 1968. D. Frydel, “One-dimensional coulomb system in a sticky wallconfinement: Exact results,”
Phys. Rev. E , vol. 100, p. 042113,Oct 2019. M. A. Miller and D. Frenkel, “Competition of percolation andphase separation in a fluid of adhesive hard spheres,”
Phys. Rev.Lett. , vol. 90, p. 135702, Apr 2003. G. Foffi, C. D. Michele, F. Sciortino, and P. Tartaglia, “Scalingof dynamics with the range of interaction in short-range attrac-tive colloids,”
Phys. Rev. Lett. , vol. 94, p. 078301, Feb 2005. J. N. Herrera and L. Blum, “Sticky electrolyte mixtures in thePercus-Yevick/mean spherical approximation,”
The Journal ofChemical Physics , vol. 94, no. 7, pp. 5077–5082, 1991. L. Blum, M. L. Rosinberg, and J. P. Badiali, “Contact theoremsfor models of the sticky electrode,”
The Journal of ChemicalPhysics , vol. 90, no. 2, pp. 1285–1286, 1989. D. A. Huckaby and L. Blum, “Exact results for the adsorption ofa dense fluid onto a triangular lattice of sticky sites,”
The Jour-nal of Chemical Physics , vol. 92, no. 4, pp. 2646–2649, 1990. J. A. Greathouse and D. A. McQuarrie, “Conventional hyper-netted chain force calculations for charged plates with adsorbingcounterions,”
Journal of Colloid and Interface Science , vol. 181,no. 1, pp. 319 – 325, 1996. J. A. Greathouse and D. A. McQuarrie
The Journal of PhysicalChemistry , vol. 100, no. 5, pp. 1847–1851, 1996. M. S. Wertheim, “Fluids with highly directional attractiveforces. IV. Equilibrium polymerization,”
Journal of StatisticalPhysics , vol. 42, pp. 477–492, feb 1986. Y. Wang, Y. Wang, D. R. Breed, V. N. Manoharan, L. Feng,A. D. Hollingsworth, M. Weck, and D. J. Pine, “Colloids withvalence and specific directional bonding,”
Nature , vol. 491, pp. 51–55, nov 2012. M. N. Tamashiro, Y. Levin, and M. C. Barbosa, “Donnan equi-librium and the osmotic pressure of charged colloidal lattices,”
The European Physical Journal B-Condensed Matter and Com-plex Systems , vol. 1, no. 3, pp. 337–343, 1998. N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H.Teller, and E. Teller, “Equation of state calculations byfast computing machines,”
The Journal of Chemical Physics ,vol. 21, no. 6, pp. 1087–1092, 1953. D Frenkel and B Smith,
Understanding Molecular Simulation .New York: Academic Press, 1996. T. Hill,
Statistical Mechanics: Principles and Selected Applica-tions Dover Publications Inc . New York, 1987. D. A. McQuarrie,
Statistical Mechanics . Harper’s ChemistrySeries, New York: HarperCollins Publishing, Inc., 1976. W. Ebeling, “Zur theorie der bjerrumschen ionenassoziation inelektrolyten,”
Zeitschrift f¨ur Physikalische Chemie , vol. 238,no. 1, pp. 400–402, 1968. H. Falkenhagen and W. Ebeling,
Equilibrium Properties of Ion-ized Dilute Electrolytes . Academic Press, New York, 1971. Y. S. Jho, S. A. Safran, M. In, and P. A. Pincus, “Effectof charge inhomogeneity and mobility on colloid aggregation,”
Langmuir , vol. 28, no. 22, pp. 8329–8336, 2012.
A. P. dos Santos, M. Girotto, and Y. Levin, “Simulations ofcoulomb systems confined by polarizable surfaces using periodicgreen functions,”
The Journal of Chemical Physics , vol. 147,no. 18, p. 184105, 2017.
S. H. Behrens, D. I. Christl, R. Emmerzael, P. Schurtenberger,and M. Borkovec, “Charging and aggregation properties of car-boxyl latex particles:experiments versus dlvo theory,”