{'fraction_non_alphanumeric': 0.04231018528866141, 'fraction_numerical': 0.01660669184863995, 'mean_word_length': 4.853984037681538, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Coupling between bulk- and surface chemistry in suspensions of charged colloids

The ionic composition and pair correlations in fluid phases of realistically salt-free charged colloidal sphere suspensions are calculated in the primitive model. We obtain the number densities of all ionic species in suspension, including low-molecular weight microions, and colloidal macroions with acidic surface groups, from a self-consistent solution of a coupled physicochemical set of nonlinear algebraic equations and non-mean-field liquid integral equations. Here, we study suspensions of colloidal spheres with sulfonate or silanol surface groups, suspended in demineralized water that is saturated with carbon dioxide under standard atmosphere. The only input required for our theoretical scheme are the acidic dissociation constants pKa, and effective sphere diameters of all involved ions. Our method allows for an ab initio calculation of colloidal bare and effective charges, at high numerical efficiency.


I. INTRODUCTION
Predicting the structural correlations in suspensions of charged colloidal particles without any fitting parameters still represents a formidable challenge of statistical physics. This is mainly due to two reasons: first, the Coulomb interactions are long-ranged and there are nontrivial correlations between the colloidal macroions and between the microions which require an extension of standard mean-field theories of linear screening [1][2][3]. Second, the (bare) charge of the colloidal particles in suspension is not known a priori, but underlies the chemical charge regulation process, with the dissociation degree of ionizable colloidal surface groups depending on the amount of added electrolyte ions and on the colloidal concentration. The resulting colloidal bare charge largely differs from the titration charge, i.e., the maximal possible charge for a colloidal particle with fully dissociated acidic surface groups [4,5].
In addition, presence of microions with non-mean-field like distributions in narrow diffusive layers about the colloidal particle's surfaces [6][7][8][9][10][11][12][13] causes that the effective electrostatic interaction is further reduced. For instance, in a one component macroion fluid model, where the microion degrees of freedom are integrated out [14], it is an effective colloidal charge that dictates the paircorrelations among colloidal particles. Typically, the colloidal (effective) charge is treated as a fit parameter. An example is the fitting of a Debye-Hückel potential to the far-field numerical non-linear Poisson-Boltzmann solution [15]. The so-determined type of effective charge is also known as renormalized charge. In experimental analysis, an effective charge is commonly used in describing colloidal static structure factors or radial distribution * mheinen@thphy.uni-duesseldorf.de functions [16][17][18][19][20][21][22][23][24][25]. Also the phase behavior [26] and the elastic properties in the solid state [5,27] can be interpreted in terms of a Debye-Hückel potential, based on a fitted effective charge, and, furthermore, colloidal effective charges determine the suspension's electro-kinetic properties [4,28,29]. Conductivity measurements, in particular, access the number of uncondensed, freely moving counterions [30,31]. Although the various effective charges, probed by these different experiments, are conceptually different from each other, the ratio of their numerical values seems to be correlated [29,31,32]. The chemical and experimental boundary conditions for charged sphere suspensions can be varied over a wide range [33], allowing for large variations in the colloidal charge numbers. In a self-consistent parameter-free approach, the colloidal bare and effective charges in an aqueous solvent should be predicted based on the chemical equilibrium conditions of dissociated surface ionic groups and bulk ions [34].
Nonlinear screening theories [35,36], computer simulations of the primitive model [37][38][39][40][41] and liquid integral equation theory of strongly coupled Coulomb systems [13,42,43] are routinely used to treat the ionic correlations, but the second aspect of bare charge variability has often been ignored in these approaches.
Monte Carlo [44][45][46][47][48][49] or Molecular Dynamics [50,51] computer simulations with an explicit account for charged surface groups are computationally very expensive, especially when the size-and charge disparity between macroions and microions is large. This renders the development of computationally more efficient methods desirable [52]. For a recent review on surface charge regulation in biomolecular solutions, we refer to Ref. [53].
Behrens, Borkovec and Grier have solved the problem of charge regulation of two electrolyte-immersed surfaces, with Poisson-Boltzmann microion distributions [54][55][56][57] and, recently, the conductivity of charged, electrolytefilled fluidic nanochannels has been investigated in a com- Fig. 1: Schematic diagram of the coupled physicochemical problem. Chemical association-dissociation balances in bulk suspension and at the colloidal particle's surfaces couple to the electrostatic and statistical-mechanical problems of variable colloidal bare charge and particle pair-correlations described in the liquid integral equation approach. Relations between the four subproblems, that are taken account of in the present work, are indicated by arrows. A closed graph of subproblems is obtained that can be self-consistently solved. parable mean-field-level study [58]. Coupled surface and bulk chemistry in colloidal suspensions has also been considered in a mean-field-like approach [59][60][61][62][63] which takes account of macroion correlations only within (revised versions of) the minimalistic cell model [15]. This was used to predict electrokinetic properties of aqueous suspensions. An account of water self-dissociation and carbon dioxide based contaminations yielded an improved agreement with experimental data in these studies on the mean-field level.
In this paper we tackle both problems -the non-meanfield correlations in ionic colloidal suspensions and the chemical regulation of the colloidal charge -simultaneously, in a self-consistent semi-analytical approach based on liquid integral equations. Thereby, ionic correlations beyond the linear screening theory level are incorporated in a good approximation. At the same time, the liquid integral equation solution provides a coupling between the chemical association-dissociation balances of acidic groups on the colloidal sphere's surfaces, and the bulk concentrations of all ionic species. The key idea, illustrated schematically in Fig. 1, is that liquid integral equations predict the excess chemical potentials of all ionic species, which then enter into the chemical associationdissociation balance of colloidal acidic surface groups. The degree of surface-group dissociation is directly proportional to the colloidal bare charge, which, in turn, influences the overall (bulk) ionic composition and the pair correlations among all ion species in the liquid integral equation system. The so-obtained implicit set of physicochemical equations is numerically self-consistently solved, yielding results that include colloidal bare and effective charges, and the suspension's pH-value.
Note that the major difficulty in tackling the coupled equation set lies in the numerical solution of the involved liquid integral equations. When all ion species are treated on equal footing in the so-called primitive model, as done in the present work, very large asymmetries between the (effective) hard-core diameters and charge numbers of macro-and microions must be resolved. These asymmetries pose a formidable challenge for the numerical stability and efficiency of solution methods for liquid integral equations. Solving the equations that occur in the present study within reasonable program execution times has been rendered possible only recently, with the advent of a numerical solution method by part of the present authors [13]. This method is based on earlier work by different groups [64][65][66][67][68], the key ideas of which have been generalized and combined in a versatile way.
This paper is organized as follows: In Sec. II, we explicate our theoretical scheme, including associationdissociation balances between all relevant reactive species in Sec. II A, constraints on the number concentrations in Sec. II B, ion pair-correlations and activities in Sec. II C and Appendix A, the effective charge number of colloidal spheres in Sec. II D, and the self-consistent solution of the coupled physicochemical equation set in Sec. II E and Appendix B. Results predicted by our theoretical scheme are presented in Sec. III, beginning with a discussion of macro-and microion pair-correlation functions in Sec. III A, and macroion bare and effective charges as well as the suspension's pH-value in Sec. III B. In the final two sections IV and V, we mention possible future continuations and extensions of the present work, and give our concluding remarks.

II. THEORETICAL SCHEME
In the following, we investigate aqueous suspensions of monodisperse colloidal spheres in thermodynamic equilibrium. Each colloidal sphere carries a mean (timeaveraged) electric charge of magnitude Ze, where Z is the colloidal bare charge number, and e denotes the proton elementary charge. In the model description applied here, colloidal spheres acquire their electric charge solely by the dissociation of acidic surface groups that are covalently bound to the sphere surfaces. We have limited our studies to two types of colloidal spheres, with either strongly or weakly acidic surface groups. The first type represents spheres that are covered with strongly acidic sulfonate (R-O-SO 3 H) surface groups. Such particles have been synthesized and used in various experimental studies of phase behavior, equilibrium and nonequilibrium properties [24,31,69]. The second type represents spheres covered with weakly acidic silanol (SiOH) surface groups, such as the experimentally frequently used colloidal silica particles [18,20,21,26,70]. Weakly acidic surface groups allow for a considerable variation of the colloidal charge by altering the suspension parameters. Both kinds of surface groups are monovalent acids.
As a consequence, we have 0 ≥ Z ≥ −N .
We denote the number concentration of particles of species i by [i], and all number concentrations in this paper are given in units of M = 1 mol/liter. The (bulk) number concentration [i] is defined as the total number of particles of species i, divided by the total system volume. In the primitive model (PM) description applied here, spherical colloidal macroions as well as monovalently charged microions are approximated as non-overlapping hard spheres with pairwise additive hard-core diameters. The methods presented in this paper could in principle be applied to suspensions including multivalent lowmolecular-weight microions, where effective charge inversion of colloidal spheres has been observed [71][72][73][74][75][76][77][78][79][80]. However, for the sake of simplicity we limit ourselves here to suspensions with monovalent microions. Throughout our analysis, we assume the approximate microion effective sphere diameters σ H3O + = σ OH − = 0.9 nm and (1) reminiscent of ions dressed with one hydration layer of H 2 O molecules. All results presented here are for colloidal spheres with hard core diameter A. Association-dissociation balances In Fig. 2, the chemical formulas of the seven reactive species of interest are given, including water (H 2 O), carbon dioxide (CO 2 ), bicarbonate (HCO − 3 ), hydronium (H3O + ), hydroxide (OH − ), and colloidal surface groups with (SgH) or without (Sg − ) an attached proton. For the systems studied in the following, SgH either stands for one sulfonate (R-O-SO 3 H) or one silanol (SiOH) group. The species in Fig. 2 are grouped by three ellipses, each surrounding the reactants of one of the three fundamental association-dissociation balances of the system, which are Here, pK SgH a and pK CO2 a are the acid dissociation constants of the surface groups and of carbon dioxide, respectively, and pK is the water self-dissociation constant. Note that Eq. (4) is short-hand notation for the combined two reactions CO 2 + H 2 O ⇌ H 2 CO 3 and H 2 CO 3 proceeding via the intermediate species carbonic acid (H 2 CO 3 ). Carbonic acid molecules are electrically neutral, and therefore do not influence the PM ion pair-correlation functions, discussed further down in subsection II C. Also, H 2 CO 3 molecules do not directly participate in either of the two reactions in Eqs. (5) and (6). It is thus unnecessary to include carbonic acid molecules explicitly into our description.
The CO 2 dissociation constant pK CO2 a quantifies the equilibrium thermodynamic activity ratio [81,82] for the reaction in Eq. (4). In Eq. (7) and further down this text, the thermodynamic activity, a i , of species i is defined according to the the convention where β = 1/(k B T ) with Boltzmann constant k B and absolute Temperature T , and where µ i denotes the chemical potential of species i. The latter can be written as the sum of the excess chemical potential µ exc i = ln(γ i )/β and the ideal chemical potential is the activity coefficient, and Λ i is the thermal de Broglie wavelength, which is of no relevance in the following. Employing the conventions in Eq. (8) implies that the reference state of substance i is an ideal gas at number density [i], with chemical potential µ id i . In the following, we approximate γ i = 1 for all electrically neutral species i. Then, Eq. (7) can be re-written as and the analogous equation quantifies the equilibrium state of the water self dissociation reaction in Eq. (5), with K = 10 −pK M. The equilibrium state of the acidic surface group dissociation reaction in Eq. (6) is characterized by Realizing that pH = − log 10 (a H3O + ), , Eq. (11) can be converted into the Henderson-Hasselbalch equation quantifying the chemical regulation of Z. Links between the three Eqs. (9), (10) and (12) are provided by the hydronium ion concentration, [H 3 O + ], and also by the four activity coefficients γ H3O + , γ OH − , γ HCO − 3 , and γ Sg − , each of which depends on the charge and concentration of all ionic species in suspension. In the self-consistent PM solution scheme used here, the intricate relations between the γ i 's, [i]'s and Z are resolved within the HNC approximation (c.f., Sec. II C).
All results presented in the following have been obtained using the (acid) dissociation constants pK CO2 a = 6.5, pK SiOH a = 4.0, and for water, carbon dioxide, silanol, and sulfonate dissociation, respectively. Values in Eqs. (15) and (16) were chosen as typical representative cases of a weak and a strong acid.

B. Concentration constraints
Without further constraints, the three Eqs. (9), (10) and (12), containing five number concentrations, four different activity coefficients, and the unknown charge number Z, do not possess an unambiguous solution. In the following, we construct a closed set of equations with a unique solution by identifying the relevant concentration constraints, and by providing PM-HNC expressions for the activity coefficients. We begin by identifying the known and unknown quantities, listed in Tab. I.
In the left column of Tab. I, the relevant known input parameters are listed, beginning with the number, N , of dissociable surface groups per colloidal sphere. This quantity is assumed to be known since, in typical experiments, it can be accurately determined by titration [30,83,84]. Likewise, the number concentration of colloidal spheres, [Col], is assumed to be known since it is Tab. I: The known input parameters of our suspension model, and the basic quantities that need to be determined for an unambiguous solution.
[Col] is the number concentration of colloidal spheres, each carrying N surface groups, |Z| of which are dissociated in equilibrium. Further explanations are given in subsection II B

Known:
Unknown: an experimentally rather well-controlled quantity. It can either be measured directly [5,26,85], or it can be calculated, e.g., on basis of a colloidal form-factor measurement, the colloidal sphere mass density, and the colloidal mass fraction [24]. In presenting our results for different values of [Col] in Sec. III, we use the colloidal volume fraction as a control parameter, since φ is more intuitively interpreted than the quantity [Col]. In Eq. (17), φ is the fraction of the total suspension volume that is occupied by colloidal spheres. Since water molecules are the overwhelming majority species, it is a good approximation to assume a constant [H 2 O] = 54.2 M, which corresponds to the number concentration of pure water. This concentration is many orders of magnitude higher than that of any other species in the self-consistent solutions reported in Sec. III.
As regards carbon dioxide, we assume a concentration of [CO 2 ]= 1.52 × 10 −5 M, which corresponds to CO 2saturated, salt-free water under an atmosphere with a CO 2 partial pressure of 3.9 × 10 −4 atm. [86,87]. Note here again that γ CO2 is equal to one in our approximate description. It is therefore consistent to prescribe the number concentration of CO 2 .
In addition to fixing [H 2 O] and [CO 2 ], a constraint arises from requiring global electroneutrality of the suspension, which can be written as The global electroneutrality constraint in Eq. (18), combined with Eqs. (9) and (10), gives the quadratic equation At this point we have collected the four Eqs. (9), (10), (12) and (19). In combination with the HNC scheme so-lution, from which the ion activity coefficients are obtained, these equations are sufficient to determine all eight unknowns listed in the right column of Tab. I.

C. HNC scheme
We employ the liquid integral equation formalism to compute the pair-correlations among all ionic species in suspension, based on the multicomponent Ornstein-Zernike (OZ) equations [88] (20) which are valid for a homogeneous and isotropic, threedimensional fluid mixture. In Eq. (20), the c i,j (r) and h i,j (r) = g i,j (r) − 1 are the partial direct and total correlation functions, respectively, between ions of species i and j.
Here, we solve the coupled OZ equations for a system of five ionic species: Number one to four are the species H 3 O + , HCO − 3 , OH − and Col, the latter denoting entire colloidal spheres that carry a charge of Ze each. The fifth ionic species is identified by the lower index 'dilCol' in the following, and represents an ultradilute fluid of colloidal spheres with diameter σ Col and with a charge of (Z − 1)e. Species dilCol is introduced merely as a bookkeeping device, necessary for the determination of the surface group excess chemical potential µ exc Sg − , as explicated in Appendix A. The number concentration [dil-Col] is selected several orders of magnitude smaller than [Col]. Hence, species dilCol exerts a negligible influence on the mutual pair-correlation functions between the four species H 3 O + , HCO − 3 , OH − and Col. To obtain a closed set of integral equations, Eqs. (20) are combined with the approximate HNC closure relation [88,89] in which the u i,j (r) are the dimensionless pair-potentials of direct interaction between ions, invoking the solvent-characteristic Bjerrum length L B = e 2 /(ǫk B T ) in Gaussian units and the pairwise additive hard core diameters σ i,j = (σ i +σ j )/2. In all calculations with results presented here, we have used L B = 0.701 nm, corresponding to water at room temperature. Assuming pair potentials of the kind of Eq. (22) for the microion and macroion species, amounts to an approximate treatment of the ion pair-interactions within the PM. The PM description neglects short-ranged van der Waals attraction, as well as changes in water polarizability which can play a role a high surface potential [90]. Furthermore, it is assumed that the charge of a colloidal sphere is homogeneously smeared out over the sphere surface. Our model thus neglects all effects arising from charge patchiness [50], a topic that has recently received much interest in studies based on the nonlinear and anisotropic Poisson-Boltzmann equation [91][92][93]. Surface charge patchiness could in principle be included into our description, if the OZ Eqs. (20) were replaced by a reference interaction site model [88,[94][95][96] description, or by anisotropic OZ equations [97]. However, the strong charge-and diameter asymmetry between macroions and microions renders already the solution of Eqs. (20)- (22) into a tedious task [13,98].
We solve Eqs. (20)-(22) by means of our recently developed method [13], which is specially well-suited for application to highly asymmetric electrolytes, in an arbitrary number of spatial dimensions. For details of the solution method, which relies on a generalized version of Ng's fixed point iteration scheme [64] and a Fourier-Bessel transform on computational grids with logarithmic spacing [65,67,68], we refer to our comprehensive description in Ref. [13]. Note here that essentially the same numerical method has been used already in the year 1980 by Rossky and Friedman [66]. Our algorithm constitutes an optimization and generalization of this earlier work, and the first application of the method to highly asymmetric electrolytes.
Once that Eqs. (20)- (22) have been solved for a given set of [k]'s and a given Z, the correlation functions are used as input for computing the thermodynamic activity coefficients of all ionic species by means of the Hansen-Vieillefosse-Belloni equation [88,[99][100][101][102][103] ln(γ i ) = βµ exc The surface group excess chemical potential, ln(γ Sg − ), which is the essential quantity in colloidal surface charge regulation described by Eq. (12), is obtained within the PM as the right-hand-side of Eq. (A5). It is taken as the sum of the colloidal sphere Coulomb self-energy change, caused by the dissociation of one surface group, plus the difference between the excess chemical potentials of colloidal spheres with charges (Z − 1)e and Ze. A brief discussion is in place here, regarding the accuracy of Eq. (23), which is the HNC approximation of an exact expression that has been derived by Kjellander and Sarman [104] and Lee [105] (see also Ref. [106]). In Refs. [105] and [107] it has been shown and discussed that Eq. (23) generally provides a very poor approximation for the excess chemical potential of particles with a hard core. Since we are indeed concerned with particles that exhibit hard-core plus Coulomb interactions, the applicability of Eq. (23) may therefore be questioned. However, our method for calculating the salient surface group excess chemical potential is based on the difference µ exc dilCol − µ exc Col between excess chemical potentials of colloidal spheres that differ in their electric charges, but not in their hard core diameters. As we have checked, the inaccurate hard-core contributions (i.e., the contributions to the integrals in Eq. (23) for 0 < |r| < σ Col, HCO − 3 ) are practically identical for both species Col and dilCol, and therefore cancel out nearly perfectly when the difference is taken. The remaining non-overlap parts of the integrals in Eq. (23) are quite accurate due to the very rapid decay of the (neglected) bridge function at non-overlap distances of particles with Coulomb interactions.
In addition to the surface group excess chemical potential, the hydronium ion excess chemical potential ln(γ H3O + ) enters into the charge regulation Eq. In a future extension of the present work, the HNC closure may be replaced by a thermodynamically partially consistent closure relation. Here, a specially suitable candidate is the closure that has been proposed by Bomont and Bretonnet [108], and that has been supplemented by an expression for the excess chemical potential [107,109], similar in form to Eq. (23), but significantly less suffering from an inaccurate hard-core contribution. Bomont and Bretonnet's closure is especially well suited for application to a restricted PM of electrolytes containing microions only, or for electrolytes containing rather small polyions like, e.g., charged globular proteins [23]. Note, however, that the application of a thermodynamically self-consistent closure to a PM with strong charge-and size-asymmetries is somewhat hampered by the fact that the number of correlation functions raises more quickly than the number of consistency criteria when the number of species is increased [110]. Therefore, keeping in mind the slight inaccuracy of Eq. (23), we resort to the simpler HNC scheme in the present work.

D. Colloidal effective charge
In the analysis of experiment results, and in the construction of theoretical schemes for colloidal dynamics, one is often interested in a mesoscopic description of reduced complexity, where the microion's degrees of freedom have been integrated out. In such a one-component macroion fluid (OMF) description, the colloidal spheres remain as the only species whose correlations are explicitly resolved, and the hard-sphere Coulomb pair-potential among macroions, u Col, Col (r), must be replaced by an effective, state-dependent macroion pair potential u eff Col, Col (r) that takes implicit account of the presence of microions.
Having solved the coupled PM-HNC Eqs. (20)- (22) for all ionic species, an effective macroion pair potential can be extracted via an inversion of the HNC relation [13,111]. In a very similar way, HNC inversion has been used to extract effective macroion potentials from digital video microscopy data [112]. The effective macroion potential from HNC inversion can be mapped to the electrostatic repulsive part, of the Derjaguin-Landau-Verwey-Overbeek (DLVO) pair potential between two finite-sized macroions in an electrolyte with microion correlations treated in Debye-Hückel approximation [14]. For the chemical composition of suspensions studied here, the square of the inverse exponential screening length κ in Eq. (24) is given by In case of a dilute suspension of weakly charged macroions with |L B Z/σ Col, Col | ≪ 1, the potential in Eq. (24) accurately represents the effective macroion pair potential with Z eff = Z. In suspensions where |L B Z/σ Col, Col | 1, the potential in Eq. (24) remains to be a good approximation of u eff Col, Col (r) at sufficiently large macroion separation distances, but the effective charge number, Z eff , satisfying |Z eff | ≤ |Z|, can considerably differ from the bare charge Z [6-12, 15, 113-123].
We determine Z eff in the following by fitting u DLVO Col, Col (r) to u eff Col, Col (r) at large particle separations. Here, Z eff is used as the only tunable fit parameter. The effective charge Z eff can be regarded as the overall charge of a colloidal sphere and that part of it's surrounding double layer in which the Debye-Hückel approximation of microion distributions breaks down.
We note here that our description of the electric double layer is similar, but not equal to the so-called 'Basic Stern Model' or 'Zeroth-order Stern Model' [124,125]. Like these variants of the Stern model, our PM description takes account of the finite size of microions in using the pairwise additive ion hard-core diameters σ i,j . However, going beyond the Stern model, our description also takes account of non-mean-field (PM-HNC) correlations between all ion species, regardless of the ion separation distance. The Stern model, in contrast, assumes mean-field (Poisson-Boltzmann) microion distributions in the diffusive (non-condensed) part of the double layer, in the same fashion as the historically preceding Gouy-Chapman model.

E. Self-consistent solution
We solve the set of Eqs. (9), (10), (12) and (19)-(23) for the eight unknown quantities in Tab. I, by the iterative algorithm described in Appendix B. This algorithm seeks a fixed point solution of the coupled set of equations by stepping repeatedly through the loop of subproblems that is schematically depicted in Fig. 1.

A. Ion pair-correlations and pH value
As a first result, Fig. 3 features the PM-HNC solutions for the partial rdf's g i,j (r) between the four ion species Col, H 3 O + , HCO − 3 and OH − , in two different colloidal suspensions, corresponding to the two panels of the figure. The partial rdf's between the ultradilute colloidal sphere species 'dilCol' and other species are indistinguishable from the corresponding functions for species 'Col', on the scale of Fig. 3, and are therefore not shown. Results in the top panel of Fig. 3 are for a suspension of colloidal spheres that carry N = 5000 silanol surface groups. In the self-consistent solution of the physicochemical problem, only 4% of the silanol surface groups are dissociated under these conditions, which results in a colloidal bare charge of Z = −212. Due to the strong electrostatic repulsion and the relatively high colloidal volume fraction of φ = 5%, the pair correlations between colloidal spheres in this suspension are rather strong, as characterized by a macroion-macroion rdf principal maximum of g Col, Col (r ≈ 6a Col ) = 2.06 (black solid curve in the upper panel of Fig. 3). The concentration of positive hydronium ions close to the negatively charged colloidal sphere's surfaces is 3.3 times higher than the suspensionaveraged hydronium ion concentration, as indicated by the contact value g Col, H 3 O + (σ Col, H 3 O + ) ≈ 3.3 (red solid curve in the top panel of Fig. 3).
The lower panel of Fig. 3 features the partial rdf's for a dilute suspension, at a colloidal volume fraction of φ = 10 −5 . Here, each colloidal sphere carries N = 230 sulfonate surface groups. Due to the small value of the surface group acidic dissociation constant, pK R-O-SO3H a = 1.5, the self-consistent solution of the physicochemical set of equations predicts that 98% of the sulfonate groups are dissociated here, resulting in a colloidal bare charge of Z = −225. Attraction of diffusing hydronium counterions towards the colloidal sphere's surfaces is strong, as signaled by the contact value, g Col, H 3 O + (σ Col, H 3 O + ) = 11.4, of the macroion-counterion rdf (red solid curve in the lower panel of Fig. 3). Noting that the low-density (mean-field) approximation g Col, 7 predicts a contact value that is two times too large, we conclude that the non-mean-field character of microion distributions is  a strong effect that must not be neglected under these conditions.
In Fig.4, we display the pH-values of different colloidal suspensions, as functions of the number, N , of acidic surface groups per colloidal sphere. Red solid curves are for colloidal spheres with silanol surface groups, and black dashed curves are for colloidal spheres that carry the more strongly acidic sulfonate surface groups. Results for three different colloidal volume fractions, φ = 10 −5 , 0.01, and 0.05 are shown in Fig. 4. At the lowest volume fraction, φ = 10 −5 , the pH-value is practically independent of N . The reason is, that the amount of hydronium ions  which are released by the colloidal spheres into suspension is negligible, compared to the number of hydronium ions created in bulk suspension in the two reactions in Eqs. (4) and (5). The resulting value pH = 5.65 is a reasonable value for demineralized water that is saturated with CO 2 under standard atmosphere. As the colloidal volume fraction is increased to φ = 0.01 and 0.05, surfacereleased hydronium ions lead to appreciable drops in the pH-value. For colloids with sulfonate surface groups, the pH-value drops more rapidly (as a function of N or φ) than in case of the weakly acidic silanol groups.

B. Colloidal bare and effective charges
In Fig. 5, we plot the absolute values of Z (black thick curves) and Z eff (red thin curves), as functions of the surface group number, N . Once again, the three volume fractions φ = 10 −5 , 0.01 and 0.05 are considered. Solid curves in Fig. 5 are for φ = 10 −5 , dotted curves are for φ = 0.01, and dashed-dotted curves are for φ = 0.05. The six rightmost curves in Fig. 5 (grouped by a blue ellipse) represent results for colloidal spheres with silanol surface groups. The six curves on the left side, corresponding to sulfonate surface groups, are nearly overlapping on the logarithmic-linear scale of the main panel. The end regions of these curves at N 230 are magnified in the inset, on a linear-linear scale. As the number, N , of surface groups increases, the colloidal bare and effective charges also increase monotonically. In case of sulfonate surface groups, Z(N ) and Z eff (N ) rise more quickly than in case of silanol surface groups. Nearly all sulfonate groups are dissociated for all probed suspension parameters, resulting in |Z| ≈ N . Dissociation of the more weakly acidic  The ratio Z eff /Z of colloidal effective and bare charge decreases as a function of N , which is due to an increasing number of microions with non-Debye-Hückel like distributions. In the suspension with φ = 0.01 and N = 2750 silanol groups per colloidal sphere, only 8% of the surface groups are dissociated in equilibrium, resulting in Z = −224, and nonlinear screening leads to a further diminished value of the effective charge of Z eff = −199.
The value of |Z| decreases monotonically when φ is raised. This is due to two reasons: First, the pH-value, entering the charge regulation Eq. (12), drops with increasing φ, starting from its CO 2 -buffer controlled limit 5.65 (c.f., Fig. 4 and Ref. [47]). Secondly, increasing φ causes increasing number densities of microions that interact electrostatically with the acidic surface groups, thereby increasing the excess chemical potential ln(γ Sg − ). Both of these contributions are generally important, reaching similar magnitudes for the φ = 0.05 silanol surface group system at high values of N .
Note that for a fixed value of Z, Fig. 5 exposes a nonmonotonic dependence of Z eff on φ: In case of silanol surface groups and |Z| = 205, for example, we find |Z eff | = 198 at φ = 10 −5 , |Z eff | = 192 at φ = 0.01, and |Z eff | = 200 at φ = 0.05, i.e., an initially decreasing |Z eff |(φ) which then increases. The same effect is also observed in case of sulfonate surface groups (see here the inset of Fig. 5). In fact, also mean-field effective charge calculations show such behavior [83,119]. The observed nonmonotonicity in |Z eff |(φ) can be understood as follows: In the infinite dilution limit, the entropic gain for counterions diffusing in the bulk beats the gain in electrostatic binding energy near the colloidal surfaces. Hence, all counterions diffuse away from the colloidal sphere surfaces, and |Z eff | is (nearly) equal to |Z| for φ → 0. When φ is increased, the expected non-Debye-Hückel like distribution of microions about the colloidal surfaces sets in, resulting in a decrease of |Z eff |. When φ is further increased, global electroneutrality demands that the microion number densities in bulk solvent (i.e., far away from the colloidal sphere's surfaces) continue to increase, and the result can be a reducing electrostatic energy penalty for a counterion that diffuses from a colloidal surface into the bulk. In the bulk, the counterion itself experiences now an appreciable screening of its electric field, caused by the presence of the many other microions. A counterion with a very strongly screened electric field will ultimately behave like an uncharged hard sphere and will not condense onto the colloidal surface at all. Therefore, at high φ, |Z eff |/|Z| can rise again. C.f., here, the similar effect that has been found in simulations of protein solution at high salinity [126].
We finally note, that the present approach is similar in spirit to the determination of effective charges from elasticity experiments [127]. There the shear modulus of a randomly oriented polycrystalline colloidal solid is determined and interpreted in terms of an effective DLVO pair potential [c.f., Eq. (24)], with Z eff as the only free fit parameter. This implies an account for nonlinear screening, but furthermore also for the so-called macroion shielding effect [128], i.e., the screening of the macroion-macroion pair potential due to the presence of other macroions. Consequently, the effective elasticity charge is lower than any electro-kinetic charge measured on the very same suspension [5,27]. Within a mean-field level description, macro-ion shielding is a many body effect [129], which considerably complicates the search for suitable pairinteractions [119,120]. It becomes most important, when the range of the repulsion exceeds the nearest neighbor distance, i.e. close to the fluid-solid phase transition. It appears to vanish at strong screening, or at elevated volume fractions [130]. If macroion shielding effects are subsumed under the elasticity effective charge, the latter can be used to predict, e.g., the fluid-solid phase boundary for this suspension employing the results of Monte Carlo simulations for charged spheres interacting via a Yukawa-type pair potential [26,131,132]. Also in the present approach all electrostatic interactions are accounted for within the PM, which naturally includes the macroion shielding effect. Our effective charge number Z eff , obtained from mapping the macroion-macroion effective interaction potential to a DLVO-type pair potential, should therefore yield a suitable input for calculations of the suspension's fluid structure on the OMF level, and allow predictions for experimentally measur-able structure factors.

IV. OUTLOOK
The theoretical scheme presented here can be rather straightforwardly generalized to aqueous colloidal suspensions with added salt or other kinds of reactive electrolytes. To this end, the salinity-dependent bulk carbon dioxide concentration can be used [86,87].
Inclusion of sodium hydroxide [83] or pyridine [133,134] into the theoretical description would be particularly interesting, since it has been reported that suspensions of colloidal silica spheres exhibit a phase diagram with reentrant fluid-solid-fluid phase sequences, when either the concentration of added base or the concentration of colloidal spheres is increased [26,135]. Constructing a closed set of equations and obtaining PM-HNC solutions for all ionic rdf's in a realistic model for a colloidal suspension with added base will be somewhat more complicated than for the solvent model discussed in the present paper, due to the larger number of neutral and ionic species that will have to be taken account of.
While we have concentrated on aqueous suspensions in this work, (variations of) the presented formalism should also be applicable to the prominent problem of charge regulation in non-aqueous colloidal suspensions [18,19,136,137], which can also exhibit unusual phase sequences like crystal-fluid-crystal [19]. In non-aqueous suspensions the Bjerrum length is one to two orders of magnitude longer than in aqueous suspensions, which results in much stronger electrostatic interactions. As a consequence, tight Bjerrum-pairing of microions occurs [138,139], and nontrivial ion correlations are of great importance in the screening of colloidal sphere charges. Incorporation of non-mean-field like ion distributions in a semi-analytical theoretical framework like the present one would therefore be desirable in case of non-aqueous media. Note, however, that in non-aqueous media the mechanisms of colloidal (chemical) charge regulation are far more complex than the simple dissociation of surface groups discussed in our present work. Charging of colloidal spheres in non-aqueous media can arise from an intricate interplay of preferential surfactant adsorption, micelle formation, and dissociation of counterions from the colloidal surfaces into the hydrophilic core of micelles [140]. One future extension of the present work should be concerned with the inclusion of these charging mechanisms into the physicochemical problem set.

V. CONCLUSIONS
We have demonstrated that a set of chemical association-dissociation balances in the colloidal bulk phase and at the surfaces of colloidal spheres can be coupled by means of liquid integral equations, and that the resulting set of physicochemical equations can be effi-ciently numerically solved. The theoretical scheme introduced here allows for an ab initio calculation of colloidal bare charges Z and effective charges Z eff for fluid colloidal suspensions in a wide range of suspension parameters. As input to the theoretical scheme one needs to know only the acidic dissociation constants pK a of the involved chemically reactive species, the (effective) sphere diameters of the macroions and of all microions, the colloidal volume fraction, and the number, N , of dissociable acidic surface groups per sphere. Different from Z and Z eff , values for N can be directly and straightforwardly obtained in titration experiments and are therefore experimentally more easily accessible.
The large macroion to microion size-and charge asymmetries in typical colloidal suspensions cause a huge numerical burden in any relevant computer simulation of the primitive model. In contrast to this, the selfconsistent numerical solution of the scheme presented here takes only few minutes or less on an inexpensive personal computer, for a given set of suspension parameters. Our method is therefore well-suited for planning and analyzing experiments with charged colloidal suspensions, and to calculate primitive model pair-correlation input for theories of colloidal dynamics including electrophoresis, colloidal diffusion and rheology.
In the PM, colloidal particles are approximated as dielectric hard spheres with solvent dielectric constant ǫ, and the electric charge is assumed to be homogeneously smeared out on the colloidal sphere's surfaces. Within this model, which neglects surface-charge patchiness, a monovalent charged surface group represents nothing else than a single elementary charge that is smeared out about the surface of it's associated colloidal sphere. This allows us to construct an approximate method to determine the charged surface group excess chemical potential in consistence with the already made PM assumptions. The three-step method consists of a colloidal sphere extraction step, a charging step, and a colloidal sphere reinsertion step, as described in the following. In order to keep the suspension globally electroneutral at all steps, hydronium (H 3 O + ) counterions are taken into account.
Step 1 (colloidal sphere extraction): From the five-component PM ionic suspension described in Sec. II C, one colloidal sphere with charge Ze is extracted and placed into pure solvent, i.e., into an infinite, otherwise particle-free, dielectric continuum with dielectric constant ǫ. To restore charge neutrality of the suspension, a number of |Z| hydronium counterions are also extracted from the suspension into pure solvent (and into infinite mutual distance). The change in Gibbs free energy in step one is thus (A1) Step 2 (charging of the sphere): Inside pure solvent, one elementary charge is removed from the colloidal sphere and placed into infinite distance from the sphere. Then, the removed charge is compressed to the hydronium ion diameter σ H3O + . The change in Gibbs free energy in this step is equal to the change in Coulomb (self-)energy of the electric charge density: Step two leaves us with a colloidal sphere of charge (Z − 1)e and |Z − 1| hydronium ions in pure solvent.
Step 3 (colloidal sphere re-insertion): Insert the colloidal sphere of charge (Z − 1)e and the |Z − 1| hydronium ions from pure solvent into the fivecomponent PM suspension. The change in Gibbs free energy in this step is: where the index 'dilCol' stands for the ultradilute species of colloidal spheres with charge (Z − 1)e each. Note that, in the thermodynamic limit, none of the five ion number densities in the suspension is changed when steps 1−3 are applied. Therefore, the hydronium ion chemical potentials in step 1 and 3 are exactly equal, and we gain the expression for the total change in normalized Gibbs free energy. The second and third row in Eq. (A4) account for the insertion of a charged surface group and a hydronium ion, respectively. Considering the excess part of all quantities in Eq. (A4), we thus arrive at the expression ln(γ Sg − ) = βµ exc Sg − = βµ exc dilCol − βµ exc Col + for the charged surface group activity coefficient γ Sg − , which is required as input to the Henderson-Hasselbalch Eq. (12) for the colloidal surface charge.