Testing the relevance of effective interaction potentials between highly-charged colloids in suspension

Combining cell and Jellium model mean-field approaches, Monte Carlo together with integral equation techniques, and finally more demanding many-colloid mean-field computations, we investigate the thermodynamic behaviour, pressure and compressibility of highly-charged colloidal dispersions, and at a more microscopic level, the force distribution acting on the colloids. The Kirkwood–Buff identity provides a useful probe to challenge the self-consistency of an approximate effective screened Coulomb (Yukawa) potential between colloids. Two effective parameter models are put to the test: cell against renormalized Jellium models.


Introduction
Colloidal suspensions contain mesoscopically large particles from the nm to µm size regime, the colloids, together with small solvent and solute molecules. These molecules are orders of magnitude smaller than the colloids, and still they heavily influence the interactions between the colloidal particles. Such microscopic species also significantly affect the thermodynamics of the system. However, they often cannot be directly visualized and most experimental techniques such as small angle x-ray or neutron scattering, probe the colloidal degrees of freedom only. This stems from the wide separation of characteristic time and length scales between microscopic and mesoscopic degrees of freedom, which therefore suggests to develop effective approaches for the colloidal particles, integrating over their microscopic counterparts, that follow adiabatically.
The focus of this paper will be on the behaviour of charge-stabilized colloidal suspensions [1]- [3]. In the subsequent analysis, we will further restrict ourselves to a mean field treatment where the correlations between microions are discarded. The resulting Poisson-Boltzmann (PB) description is valid in aqueous solvents under usual conditions of temperature and for monovalent micro-ions, because existing colloids do not allow us to reach the high Coulombic couplings that are required to observe deviations from mean-field [3]. In the following, we therefore focus on the more special case of charged colloids in a simple 1 : 1 electrolyte, adopt the PB approach, and treat the micro-ions as point-like particles dissolved in a continuum of solvent molecules.
Even within such a simplified framework, the solution corresponding to N c interacting colloids is a difficult problem from a numerical point of view [4,5]. In principle, one has to compute the density distribution of microions in space ρ micro (r) for any given colloidal configuration. First one has to solve the PB equation around a large collection of N c colloids. One can then compute the corresponding stress tensor to obtain the force felt by each colloid, 4 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT

Two models for effective charges and screening lengths
We start by outlining the two effective charge models that we here are concerned with. The idea behind using effective parameters to incorporate effects of nonlinear screening in a pair-potential based otherwise on linear theory, is explained and reviewed in [10].
Although the effective potential has an unambiguous definition, there is no rigorous operational route to construct this object and it is common belief, under scrutiny here, that a Yukawa form [1]- [3] βu eff (r) = Z 2 eff λ B exp(κ eff a) 1 + κ eff a 2 exp(−κ eff r) provides a reasonable description. Here a is the radius of the colloid (assumed to be a sphere), λ B is the Bjerrum length, kT = β −1 is the thermal energy, while Z eff and κ −1 eff are the effective charge and screening length. Such a 'DLVO'-like expression [11] would accurately reproduce the large distance interaction of two colloids immersed in a salt sea [1]- [3] but it should be kept in mind that it has to fail at short distances [1]. One should also keep in mind that even if one had the perfect effective potential for two colloids, one cannot be sure that it is also the appropriate effective pair-potential for a concentrated suspension of colloids, because in general these effective pairpotentials cannot be superposed. Many-body interactions between the colloids must be taken into account [12], as they can have a considerable impact on the structure of the suspension [13]. We also remark, that in the salt-free case, the validity of (1) is not obvious. Other limitations of the concept of effective pair-potentials are discussed elsewhere [1]- [3], [14].
In appendix A, we recapitulate the essentials of the two empiric models (Jellium and PB-cell model) used here to determine Z eff and κ −1 eff in the Yukawa pair-potential. As emphasized in the introduction, effective parameters can also be used to predict the osmotic pressure of the suspension. We define the osmotic pressure from the pressure P and the reservoir pressure P reservoir , as = P − P reservoir = P − 2c s kT, (2) where c s denotes the (monovalent) salt concentration in the reservoir. Furthermore, the osmotic isothermal compressibility χ is defined via equation (2) as where the derivative with respect to the colloid density ρ c is taken at constant salt chemical potential (i.e. constant c s ). Applying these two definitions, the osmotic pressure of the suspension can be related to the effective parameters (see appendix A) with κ 2 = 8πλ B c s . The subscript 'micro' is explained further below. The compressibility then reads

5
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT in which η = 4πρ c a 3 /3 is the colloid volume-fraction. In these expressions κ 2 eff a 2 is given either by the PB-cell expression, equation While the Jellium model is implicitly based on the assumption that the colloid-colloid pair distribution function g(r) 1, the PB cell model originally takes the opposite view and assumes a g(r) corresponding to a crystalline colloidal configuration. This latter assumption rests on the observation that the repulsively interacting colloids arrange their positions such that each colloid has a region around it which is void from other colloids and which looks rather similar for different colloids. In other words, one model assumes a regular arrangement of the colloids in the suspension and thus a rather solid-like structure, while the other model assumes no structure at all-a situation typical of fluids at low density. In that sense, the two models chosen are complementary to each other.
The effective charges as a function of the bare colloid charge Z bare go to Z bare for low bare charges, and to a saturated effective charge Z sat eff for high charges, and can therefore be roughly approximated by Except in section 5, we will mostly concentrate on the saturated regime that is most suited to describe colloidal suspensions [15] (for a discussion of the crossover behaviour between low and high bare charges see [16] for symmetric electrolytes and [17] for 1 : 2 and 2 : 1 electrolytes). By choosing in all our calculations a Z bare large enough to ensure saturation of the effective charges, we can thus get rid of the additional parameter Z bare .
We have plotted the effective saturated charges as a function of η in figure A.1 of appendix A. One can identify a regime in η where salt-ions dominate the screening, leading to effective charges that are practically independent of the volume-fraction. At other values of η however counter-ions outnumber the salt-ions; this is the regime where screening is dominated by the counter-ions. As the number of these counter-ions depends on the number of colloidal particles, the effective charges become η-dependent when the counter-ions dominate the screening. One may estimate from equation (A.2) of the appendix the volume-fraction η * where one may expect to find the crossover between both screening regimes, We remark that this threshold value η * was derived within the Jellium framework. We have no clear definition like this for the cell model and we have determined the η * for the cell model empirically (see subsection 4.2). Empirically, as will be discussed in subsection 4.2, the crossover volume fraction for the cell model is found to be much smaller, η ≈ 0.2η * . Working in the regime of saturated effective charges, we are left with in total only two independent input parameters: the volume fraction η and the salt concentration κa of the salt reservoir. Multiplying Z sat eff with λ B /a and the osmotic pressure with 4πλ B a 2 , as we have done in equations (A.2), (A.5) and (4), we avoid λ B /a as an additional independent parameter. However, the colloid-colloid pair correlation function g(r) does not exhibit the same scaling behaviour as micro so that it is important to precise which value of λ B /a has been used in the simulations. Unless otherwise specified, we have considered λ B /a = 0.01, a reasonable value for colloidal systems in an aqueous environment (smaller values are also met in experiments). As for η we have stopped our calculations at high volume fraction when the liquid started to solidify, while for κa we have worked in between two extremes, κa = 1.5 (high-salt regime) and κa = 0.0 (no salt). Figure 1 compares the pressure prediction micro of both models. The curves are based on equation (4) in combination with (A.3), (A.4) and (A.2). In the no-salt case both models have the same low-dilution behaviour, but different limiting behaviour in the presence of salt ions where we have an algebraic decrease for the Jellium, but an exponential one for the PB cell model. The agreement between Jellium and PB-cell is excellent in the no-salt case, up to volume fractions around 0.1 which is remarkable in view of the differences in the effective charges.
Before proceeding we wish to emphasize that many more effective-charge approaches can be found in literature, including density functional theory (DFT)-based schemes (e.g. [18]), MC solutions of the cell model (e.g. [19,20]), and other mean-field models with various criteria for the effective charge (see reviews in [10,14]). We here have selected both the Jellium and the PB-cell models because they are complementary to each other and because they are probably the two most practical effective-charge models; they are sufficiently accurate, well-established in literature, quite simple and rather straightforward to implement (see the appendix of [21]).

Test 1: many-body forces versus Yukawa forces
We are now in the position to perform the first test. We recall that the interactions among colloids in a suspension are mediated by the micro-ions and are therefore in origin complex manybody interactions. Integrating the micro-ions out, the force on any given colloid depends on the positions of all other colloids in the system. It is not straightforward that this force can be written as a sum of pairs only. In principle, one has to sum over pairs, triple, . . . and over all many-body configurations. To test how important the higher many-body contributions are, the true many-body forces have to be compared to the forces obtained by summing up the effective pair interactions. 6 This is done in this section. We arbitrarily selected a few typical colloidal configurations from a MC simulation of 4000 Yukawa particles and solved the nonlinear PB equation in the region between the colloidal spheres of this colloidal configuration using our multi-centred PB solver described in detail in [6]. As in [6] we integrated the stress-tensor around each colloidal sphere to obtain the force acting on each colloid. These forces F MB -which include all many-body forces acting on the individual colloid-can now be compared to the forces F YU one obtains by summing pair-forces derived from Yukawa-potentials, equation (1), with effective PB-cell or Jellium parameters. We have evaluated the ratio of the magnitudes F YU /F MB and the (cosine) angles between the forces F YU · F MB /F YU F MB for each particle and have plotted the corresponding histograms, that have been averaged over a few configurations. The results of these calculations are shown in figure 2 for the salt-free case and for the high salt case, κa = 1.5.
We see from figure 2(A) that the ratio between the forces is very close to one and has a narrow delta-like distribution in the high salt case, as we should expect, since the Yukawa 8 Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT one-component picture has already proven to be applicable at such parameters. Interestingly, in the no-salt case, the force ratio shows a broader distribution but is still quite peaked, although not centred at unity, see figure 2(B). The forces based on the Yukawa one-component picture are therefore on average larger (cell) or smaller (Jellium) than the correct many-body forces. We may then anticipate that when used in the OCM, Jellium parameters will underestimate the structure, while cell data will lead to an overestimation (see figure 7(B)). Despite these deficiencies, both models seem to reproduce the forces reasonably well, which is especially true for the angular distribution displayed in figure 2(C). An interesting result is the bump observed in figure 2(B) at low force ratios: here not any of the effective-charge models, but the Yukawa pair-force model as such fails dramatically, predicting a total force that is five times smaller in magnitude than the correct many-body force. We have investigated the local structure around the colloids feeling these forces. They are always located in regions of larger local density. Here, the local mean-distance between the colloidal particles is relatively short (compared to the suspension-wide mean distance) and the Yukawa pair potential is thus probed at rather short inter-particle distance. However, at these distances expression (1) is bound to fail, as we have indicated already in section 2. The extent of this failure is now quantified in figure 2(B).

Test 2: thermodynamic consistency and the Kirkwood-Buff relation
Having seen in the previous section that the distribution of forces felt in situ by the colloids in the mixture can in most cases be captured by a Yukawa effective pair potential, we now turn to the second test of these potentials and relate the compressibility of the suspension to its structure as computed within the OCM. The way to do so is to apply the Kirkwood-Buff relation introduced further below. micro of both the PB-cell and the Jellium model is depicted in figure 1. It provides an excellent approximation to the total pressure of salt-free suspensions: both models lead to a pressure that is in very good agreement with existing experimental data [23] and PM simulations [24], 7 see e.g. [9,25]. 8 The Kirkwood-Buff relation (10) now allows for a test of the effective potential chosen: u eff should lead to a colloid structure, embodied in a long-wavelength structure factor S(0) computed within the OCM, that is compatible with the compressibility following from micro . As the same effective parameters appear both in the Yukawa-forces leading to S(0) and in micro of equation (4) this provides a critical consistency test of our effective-charge models.
We will also report results at medium and high salt concentrations, where it does not seem possible to test the effective potential as severely as in the no-salt case, but where it is still of interest to compare the different pressures introduced in this section. A similar idea has previously been worked out by Lobaskin et al [28] who checked the consistency of thermodynamic and structural properties of an asymmetric electrolyte containing macroions with 60 elementary charges and monovalent counterions. Primitive-model results have been compared with different more approximate theories such as the cell model (within mean-field and even beyond). Finally, technical details may be found in appendices B and C.

9
Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT We start this section with an outline of the theoretical background and present the actual results of our second test in subsection 4.2.

Theoretical background
4.1.1. The osmotic pressure of a charge-stabilized colloidal suspension. Our problem is best treated in a semi-grand canonical ensemble. N c colloids are confined to a certain volume V which the small solvent molecules and additional micro-ions are free to leave; they are coupled to a reservoir fixing their chemical potentials. For the moment, we need not distinguish between the different types of small molecules, but may denote them collectively by an index s, while the colloidal degrees of freedom are referred to by the index c. The total internal energy of the system can then be written as U = U cc + U ss + U cs where the index combinations indicate over what degrees of freedom the bare interactions have to be summed. The force on particle i (i ∈ c, s) then is F i = −∇ i U and the total pressure of the colloidal suspension can be obtained using the virial where the thermal average has to be taken over the colloidal as well as solvent/solute degrees of freedom. The first term, the ideal-gas term, is just the sum of densities of all components. Within the OCM, the pressure reads where the forces are obtained from the gradient of the effective potential and the summation runs over the N c colloids only. By definition, the effective potential u eff reproduces the same pair (colloid-colloid) distribution function g(r) as that of the original mixture, assuming a pair-wise Hamiltonian, see e.g. [1]. However, there is no guarantee that the corresponding sum of effective pair forces coincides, for a given colloid i in a given colloidal configuration, with the micro-ion average of −∇ i U which provides the true force. The comparison between both forces has been made in section 3. For the sake of the discussion, we assume here that they coincide. Equation (9) then uncovers one severe deficiency of the OCM picture: P OCM only provides a contribution to the total pressure given by equation (8). In particular, in the limit of low salt the micro-ions are known to dominate the overall pressure of the suspension: that is the reason why micro of figure 1 provides such an excellent approximation to the total pressure. This implies on the other hand that P OCM is a very poor approximation for P. The density derivatives of both pressures are however intimately connected, under all conditions, as we discuss in the following section (see equation (14)).

The
Kirkwood-Buff relation. The Henderson theorem [29] guarantees that every paircorrelation function g(r) can be uniquely associated with one pair-potential. As for the system considered here, this implies that for any given colloid density, there exists in principle one, and only one, effective pair-potential leading to a g(r) within the OCM that is in perfect agreement with the correct colloid-colloid pair correlation of the full multi-component system. Unfortunately, this does not guarantee that in such a multi-component system the higher-order correlation functions are equally well reproduced (see for example [13] where pair-, but not triplet correlations of a many-component colloidal systems could be reproduced within a simple pair-interaction picture).
However, there is one important quantity one can still derive if the correct colloidcolloid g(r) is known: the osmotic isothermal compressibility χ defined in equation (3). The Kirkwood-Buff relation relates this thermodynamic quantity to the infinite wavelength limit of the colloid-colloid structure factor S(q) [30] where S(0) is related to g(r) via Remarkably, equation (10) is exact: nothing about the micro-ion-micro-ion correlations or the micro-ion-colloid correlations needs to be known, but only the colloid-colloid pair correlations are required to compute the correct osmotic compressibility of the full multi-component system. By contrast, the full compressibility of the system, defined as (8), depends on the pair correlations of all components [1]. Connecting thermodynamic to structural information, equation (10) is of central importance in the present work. In many colloidal suspensions the micro-ions determine the thermodynamics, while the macro-ionic degrees of freedom are more important for the structural properties of the system. In these cases equation (10) is also well suited to bridge the gap between the micro-ion and the macro-ion oriented viewpoints, thus providing a severe test for the quality of the effective potential. Indeed, though having a clear-cut definition, an effective potential is nevertheless a difficult object to compute and one often postulates its functional form, just as we have done in equation (1). With equation (10) we can investigate a posteriori the relevance of the underlying ad hoc assumptions made in deriving the effective potentials. Regarding the effective-charge concepts under scrutiny here, equation (10) describes a route how to connect u eff with micro . As the same effective parameters appear in both quantities, this relation also provides a consistency test of our effective-charge models.

4.1.3.
On the various contributions to the total pressure. So far we have introduced three different pressures: P from which the osmotic pressure follows, P OCM , and finally P micro , the latter quantity being for instance the pressure found in the cell or in the Jellium models, see equation (4). At this point, it seems worth discussing the non-trivial connection between these three pressure expressions, at the expense of introducing a fourth quantity: we define a 'colloidal' pressure P coll from the difference between the total pressure P and P micro P = P coll + P micro .
The pressure P micro can be considered as arising from so-called 'volume terms' in the total free energy of the system, see e.g. [31,32].

Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT
In the situation where the system is confined in a closed box of volume V (ρ c = N c /V ), it can be shown that 9 ,10 The third term on the right-hand side in (13), for which the surface integral with normal oriented outward runs over the box confining the system, accounts for the direct coupling between colloids and micro-ions. It is precisely this quantity that the micro-ionic cell-model approaches aim at computing. In other words, this third contribution may be identified with P micro in (12) so that remembering equation (9), we obtain P coll P OCM . However, a closed cell is not the most convenient configuration (in particular, the effective potential of interaction not only depends on the relative distance between two colloids, but also on the distance to the confining walls), and from a computational perspective, it is often more desirable to work with periodic boundary conditions systems. Unfortunately, equation (13) which provides a physically transparent interpretation to P coll in the confined case (close cell), breaks down with periodic boundary conditions. This failure is discussed at length in [33], together with the fact that in low salt (or no salt) conditions, one has P coll P micro which implies P P micro (see [34] for a related discussion).
From the previous discussion, it appears that the connection between P OCM and P is not straightforward. The simplest relation between both quantities follows from the remark that the colloidal structure within the OCM is of course the same as in the original mixture (assuming the effective potential to be the correct one). Equations (10) and (11) then dictate that the compressibilities in both approaches coincide: where it is crucial to compute the OCM compressibility at constant potential of interaction u eff , i.e. discarding any density dependence of the effective potential. In a region of parameter space where the density dependence of u eff is absent or weak enough, the 'salt' and 'potential' subscripts in the partial derivatives of equation (14) correspond to the same constraint, since u eff then only depends on the salt chemical potential, besides relative distance. It is then possible to integrate equation (14) to obtain P P OCM . Such a situation is met in the salt dominated regime to be discussed below.

Applying the Kirkwood-Buff relation
We next apply equation (10) to the case at hand and expect to find the following: In the nosalt case where the total pressure is very close to P micro we should obtain within the OCM a long-wavelength structure factor S(0) compatible with ∂ρ c /∂(βP micro ). Any difference between both quantities indicates a lack of consistency in the (somewhat uncontrolled) procedure leading to the effective Yukawa potential (1) chosen here. Such an inconsistency may take its roots in an incorrect computation of effective parameters (charge and screening length) or more fundamentally in the functional screened Coulomb form taken.
We have calculated S(0) within the OCM as a function of η as described in detail in the appendices B (MC simulation) and C (integral-equation theory). Both methods lead to identical results, shown in figure B.1 of appendix B. Such sets of S(0) curves had to be computed twice, using once the PB cell model parameter and once the Jellium model parameter within the effective Yukawa pair-potential. Figures 3 (PB-cell) and 4 (Jellium) show the corresponding S(0) curves. Each figure consists of three graphs, corresponding to κa = 1.5, κa = 0.5 and κa = 0.0. In each graph, we compare S(0) to the compressibility χ micro following from the ρ c derivative of pressure P micro (or equivalently micro ). To distinguish between the two effective charge models used to compute micro via equation (4), we introduce the notation χ micro (PB) and χ micro (Jell). Alternatively, one may compute the OCM pressure P OCM , defined in (9), and obtains then the compressibility χ OCM as the ρ c derivative of P OCM . If PB-cell model (Jellium model) parameters have been used in the effective Yukawa potential, we denote the resulting compressibility by χ OCM (PB) (χ OCM (Jell)). These quantities, and in addition the inverse of the saturated effective charges, are also graphed in figures 3 and 4.

No salt.
We observe that χ micro , in the no-salt case, is close to S(0). The agreement seems to be slightly better for Jellium compared to cell parameters at low density, while the opposite holds at higher densities. This finding is consistent with the fact that the probability distribution shown in figure 2(B) (where η = 10 −4 , a low value) is slightly more peaked for the Jellium than for the cell. On the other hand, the inadequacy of χ OCM for κa = 0 is expected: P OCM has here nothing to do with the total pressure P; the resulting compressibility is about a factor of two larger than S(0). This question is addressed in [33] and for completeness, we adapt the argument in appendix D. As can be seen in figures 3 and 4, the no-salt compressibility is close to the inverse effective charge. This stems from the fact that at least within the Jellium, one has βP micro = Z eff ρ c for c s = 0 and since the density dependence of Z eff is mild (logarithmic, see [9]), one has χ micro 1/Z eff . Within the cell model, βP micro = Z eff ρ c is no longer exact, but accurate enough to lead again to χ micro 1/Z eff .

Added salt.
In the first two figures in figures 3 and 4 (the ones with κa > 0) we can see two regimes. At low volume fraction the OCM compressibility χ OCM is a good approximation and the micro-ionic χ micro fails, while at high enough volume fraction the reverse is true, χ micro is good and χ OCM fails.
The derivative involved in the computation of χ OCM includes the density dependence of the effective parameters. There are however situations where these parameters are virtually independent of ρ c (see figure A.1). This is the case for example at low density or high enough salt concentration. Looking back at equation (14), one realizes that in these cases the constant potential of interaction constraint coincides with that of constant c s . Then equation (14) dictates that χ OCM must be equal to S(0). It is therefore not a surprise, nor a deep finding, that in figures 3 and 4 the S(0) coincides with χ OCM in all those situations where the effective parameters have no ρ c dependence, i.e., at high enough salt concentration and at low density.
Such an agreement notwithstanding allows us to compute the full pressure P when the effective potential is density independent. Remembering that at high salt content, the effective forces following from u eff are very close to their 'exact' counterpart (see figure 2), we expect P OCM to provide there a good approximation to P, at least at not too high densities, i.e. in the salt-ion dominated screening regime when η η * with η * from equation (7). At high η the counter-ions dominate the screening. Thus, no matter how large κa is, there is always a regime at high enough η where the salt-ions can be neglected and where the results become independent of κa. That implies that curves differing in κa must approach the same value at high enough η. That can be observed, for example, in figures A.1(A) and B.1, where all curves show the same high η behaviour. The same applies to figures 3 and 4. At high densities, all curves approach the zero-salt curves and go to 1/Z eff , which is the same for all κa in this counter-ion dominated limit. On the other hand, for η → 0 the compressibilities for κa = 0 approach the ideal-gas limit.
We have already remarked that for κa = 0, χ micro coincides with the OCM S(0). Since, for large η (i.e. η > η * ) all κa = 0 curves must approach the κa = 0 curve, we expect to find that at high enough η, χ micro coincides with S(0), since then the system is close to the salt-free regime. 11 This is indeed the case in figure 3 (and to a lesser extent in figure 4 where only the trend is visible). This limiting behaviour seems to be a natural self-consistency requirement to impose to any effective potential.
For the Jellium model, we have already defined the crossover volume fraction η * (equation (7)) separating the salt and counter-ion dominated screening. The dependence of η * on the salt concentration κa and on the volume fraction η is graphed in figure 5 (solid line). This graph is meant to summarize our findings of figures 3 and 4. For the Jellium model, the counterion dominated regime (η > η * ) is to the right of the solid line. Here, the OCM S(0) is found from figure 4 to be consistent with χ micro . The other regime (the salt-ion dominated regime) where we can observe from figure 4 that S(0) is consistent with χ OCM is unfortunately not located to the left of the solid curve, but rather to the left of the dashed line empirically determined as η * s = 0.2η * . This means that in frame of the Jellium model there is a gap between the salt dominated region and counter-ion dominated region (η * s < η < η * ). At these intermediate volume fractions neither Such a gap is not found using parameters from the PB-cell model. For the PB cell model we have no analytical prediction for the value of the crossover volume fraction and have to read the values from the results in figure 3. Empirically we have found the crossover volume fraction 11 In the salt-free case, the suspension behaves as if it were permanently in the 'large' η-limit and S(0) shows a weak η dependence only because of the η-dependence of Z sat eff . This special feature of the zero-salt case can also be appreciated from the pair-correlation functions displayed in figure A.1(B). Plotting these functions over the distance in units of the mean separation, they largely resemble each other for all η considered, indicating a degree of particle correlation that is almost independent of η.  (7), while on the dashed line η is equal to 0.2η * .
of about η ≈ 0.2η * , which is the same as the previously defined η * s for the Jellium model. This is plotted as a dashed line in figure 5. For the PB-cell model, we may then summarize figure 3 as follows. For system parameters lying on the right-hand side of the dashed curve, S(0) is consistent with χ micro (PB), for those system parameters on the left hand side S(0) is consistent with χ OCM (PB).
We emphasize at this point that with salt at high η, a better self-consistency is not necessarily synonymous with the fact that P micro is in itself a better approximation for P. The situation is different for κa = 0 where we have the extra piece of knowledge that P P micro , that comes from comparison with experiments or PM computations (see also section 5).
After this discussion, we can summarize the results of our Kirkwood-Buff consistency check of the effective Yukawa pair-potentials and the two effective charge models, as follows: 1. For κa = 0, we find consistency essentially for both effective-charge models where the agreement is better for the Jellium model at low densities and for the PB-cell model at high densities. 2. For κa = 0, the compressibility derived from P OCM is inadequate and not compatible with the OCM S(0). The OCM approach with the Yukawa potential as effective pair-potential seems to produce correct results for the structure of the suspension, but not for the pressure. 3. The no-salt compressibility is close to the inverse effective charge. This offers a convenient way to estimate compressibilities χ micro at κa = 0. 4. For κa = 0 and low volume fraction (η < η * s ), consistency is found between χ OCM and S(0), but not between χ micro and S(0). The OCM approach makes sense for calculating both the structure and the pressure. 5. For κa = 0 and high volume fraction (η > η * for Jellium, η > η * s for PB-cell), we get back to the zero-salt case where S(0) is consistent with χ micro , but not with χ OCM . 6. The value of η * s = 0.2η * has been empirically determined at the constant value of λ B /a = 0.01. This result might depend on the value of λ B /a, which has not been studied in scope of this paper.
In cases where we have reason to assume that the OCM S(0) is close to the correct χ, we can of course estimate the pressure P S(0) from integrating the OCM S(0) To illustrate this idea, we have performed this integration in figure 6 for κa = 0 and κa = 1.5, using the PB-cell effective parameters. The resulting curves can then be compared to micro from equation (4) (again with the PB-cell parameters), which for κa = 0 we know to be an excellent approximation to the full pressure of the suspension. We have data for S(0) starting at volume fraction η 0 ≈ 10 −8 , which is used as the lower bound of integration in equation (15). A small constant is then added to the right-hand side of (15), so that the result coincides with micro at η 0 . The figure demonstrates that P S(0) agrees well with micro , surprisingly, not only for κa = 0 but also for κa = 1.5. However, at very low η, the agreement is excellent for κa = 0, but not for κa = 1.5, as one would expect from figure 3. Note that a hypothetical discrepancy between P micro and P S(0) here could not be considered as a lack of self-consistency, since the contribution P coll in P = P coll + P micro can be non-negligible.

Test 3: comparison to PM data
P Linse has carried out a full PM simulation of a salt-free charge-stabilized colloidal suspension and computed the pressure P as a function of volume fraction η [24], see also [28]. The highest  [28]. From left to right, the packing fractions are η = 0.08, 0.04, 0.01, 0.0025, 0.00125. (C) Compressibilities derived from the data in (A) via equation (3) and from the colloid-colloid S(0) as indicated.
colloidal charge considered in this work does not lead to fully saturated effective charges (Z bare λ B /a 14.2). Figure 7(A) shows the osmotic pressure as a function of η from [24] and compares it to P micro (= micro if κa = 0). P micro is given by equation (4) with either PBcell or Jellium model data. For comparison, we also added the cell-model pressure one would obtain if the charges were saturated. This figure illustrates the quality of PB-cell and Jellium P micro in the de-ionized limit, which has been repeatedly emphasized in the previous analysis. Figure 7(C) shows the corresponding compressibilities, and as expected from figure 7(A), it can be seen that χ micro is very close to the correct PM compressibility. Figure 7(C) supports our findings of the previous sections: at low density, the Jellium effective potential leads to a S(0) that fares slightly better than its cell counterpart. The PB cell and Jellium S(0) respectively provide lower and upper bounds for the true compressibility, and correspondingly, respectively upper Institute of Physics ⌽ DEUTSCHE PHYSIKALISCHE GESELLSCHAFT and lower bounds for the pair distribution function g(r) (see figure 7(B)). This can be traced back to the force distribution seen in figure 2(B) where the Jellium effective potential leads to a slight underestimation of the 'exact' force distribution, while the opposite holds for cell data. Somehow, considering an appropriate mean of Jellium and cell effective parameters would improve the quality of the predictions. A step forward in this direction has recently been made in [35]. Of course, at very high Coulombic couplings, non-mean-field effects would spoil the scenario put forward here, see e.g. [36] for an illustration, where a phase separation is reported).

Summary and conclusion
The forces between colloids in a charge-stabilized colloidal suspension are commonly approximated by a sum of pair forces derived from the gradient of a screened Coulomb potential u eff with effective charge and screening length. A number of models can be found in literature on how to determine these effective parameters. We here computed them either from the PB cell model supplemented with the Alexander et al recipe [37,38], or from the renormalized Jellium model [9,39]. The effective parameters are therefore not fitted but derived from a well defined although difficult to control procedure. It was our main motivation to assess the relevance of this procedure, be it in the Jellium or in the cell case. To this end, we have performed three tests.
Test 1: we have performed in situ measurements of effective forces in typical colloidal configurations, from the solution of a multi-centre PB theory (with N c = 4000 colloids). We are not aware of similar measurements in the literature (see footnote 6). The distribution of these forces has been compared to forces obtained from summing Yukawa effective forces. Under high salt conditions, the pair-wise and many-body approaches gave very similar force distributions (the test was performed under conditions where cell and Jellium effective parameters are very close). In the no-salt regime, the agreement is less quantitative but still quite good given that the system is there always strongly correlated, with strongly overlapping double-layers. Jellium(/cell) forces slightly underestimate(/overestimate) the 'true' force. We also found situations where the Yukawa pair-force model fails dramatically, predicting a total force that is up to five times smaller in magnitude than the correct many-body force (bump observed in figure 2(B) at low force ratios).
Test 2: we have used the Kirkwood-Buff identity to check the consistency in the procedure leading to the effective Yukawa potential u eff and the two effective charge models. On the whole, the consistency of u eff is remarkable given the simplicity of the underlying procedures. Our findings in detail have already been summarized at the end of subsection 4.2. The most important results are: 1. In the salt-free case, we find consistency for both effective-charge models. The micro-ionic contributions χ micro are consistent with the OCM S(0). The Jellium potential is slightly better at low densities, but performs less than its cell counterpart at higher volume fraction. By contrast, the compressibility derived from P OCM is not compatible with the OCM S(0) and consequently the OCM virial pressure P OCM is far from the true pressure P. 2. With added salt, the P OCM provides a good approximation to the total osmotic pressure if the volume-fraction η is lower than the threshold value 0.2η * , (η * of equation (7)). This is the case essentially if u eff is density independent, i.e., at low enough colloid density and/or high enough salt content.
proposed in [40]; it is however an important ingredient to account for nonlinear screening effects. Z eff thus determined is used to calculate which may be rewritten as where η = 4πρ c a 3 /3 is the colloid volume fraction while κ 2 = 8πλ B c s is the squared inverse screening length in the reservoir.

A.2. The PB cell model
The PB cell model [21,37] rests on the observation that the repulsively interacting colloids arrange their positions such that each colloid has a region around it which is void from other colloids and which looks rather similar for different colloids. In other words, the Wigner-Seitz cells around two arbitrarily selected colloids are comparable in shape and volume. One now assumes that the total charge within each cell is exactly zero, that all cells have the same shape, and that one may approximate this shape such that it matches the symmetry of the colloid, i.e., spherical cells around spherical colloids. The cell radius R is chosen in consistency with the colloid volume fraction, and the PB equation within the cell is solved with appropriate boundary conditions at the cell edge and the colloid surface. Thus, through the finiteness of the cell plus the boundary conditions, the presence of all those colloids not inside the cell are taken into account. For a more detailed description of the cell model approximation see [8].
From the numerical solution of the PB equation, one obtains the electrostatic potential at the cell edge φ R , and can now proceed to compute the effective screening parameter if κ 2 a 2 > 0 and if κ 2 a 2 = 0 where µ 2 appears in the PB equation ∇ 2 φ = −µ 2 e −φ of the salt-free case and is determined from the electro-neutrality condition. Then, the effective charge following Alexander and collaborator's recipe [37] is given by [15,21,41] . There are some differences, but in the first approximation the curves superimpose. Inset: the ratio 1/(κ eff d nn ) of the double layer thickness (derived from the cell approximation) and the mean distance between particles as a function of the colloid volume fraction η at three different salt concentrations.
A simple approximation valid for large bare charges is given in [38]. We emphasize here that following such a route to define Z eff and κ eff , a 'natural' relation such as that embodied in equation (A.1) is lost, except in certain particular limits [21] (low density, or high density, or low charge). We also note that it has been shown recently that the cell model effective charge is accurately reproduced by a dynamical rule which defines the condensed micro-ions through a bound on their total energy [42], a criterion that may also be considered when mean-field breaks down. Figure A.1 serves to discuss and compare the two effective charge models in the (κa, η)parameter space. Figure A.1(A) compares the volume-fraction dependence of the saturated effective charges. The change from salt-ion to counter-ion dominated screening (occurring at η η * ) can be recognized from the onset of a η-dependence of the effective charges. For the no-salt case screening is always due to the colloidal counter-ions, and, indeed, both models predict a strong η-dependence of the effective charges even in the limit η → 0. At low η, a suspension of charged colloids with no extra salt will always be correlated regardless of how low a volume fraction is considered. The reason is that the thickness of the double layer then grows faster than the mean-distance between the particles, that is, the ratio of the double-layer thickness 1/κ eff and the mean distance d nn = ρ −1/3 c grows with decreasing volume fraction, as opposed to the case with external salt where this ratio decreases (see [1] and the inset of figure A.1(B)). Figure A.1(A) demonstrates that the two effective charge models agree inasmuch as the η → 0 limiting behaviour is concerned, 12 but disagree in the opposite limit with the Jellium model generally predicting an earlier change from the salt-ion to the counter-ion dominated screening and smaller effective charges. For high η, the curves for all three values of κa must ultimately converge for each model when the contribution of salt-ions to the screening ceases to be significant. A rather special feature of the Jellium model is the pronounced minimum in Z sat eff at intermediate volume fractions observed for all values of κa, something that is not present in the PB cell model for the salinities investigated, but that would be observed at lower salt.

Appendix B. Numerical procedures
Technical details for the calculation of the effective charges and screening parameters are given in [9,39]. Within the OCM picture, we have performed MC simulations, typically with 10 000 particles in a simulation box of side length L box applying periodic boundary conditions. We have carried out 5 × 10 6 MC cycles to reach thermal equilibrium and another 2 × 10 7 cycles for the averages. The pressure has been obtained from equation (9), the structure factor directly from and S(0) has then been approximated by S(q min ) where q min = 2π/L box . 13 Special care has been taken that the value converges with the system size. As an independent check we computed the static structure factor S(q) from the g(r) by Fast Fourier Transform and checked that the value at q min was the same. Much faster than MC simulations are structure calculations using the Ornstein-Zernike (OZ) equation [43], see appendix C. This integral equation has been solved using the well-tested Rogers-Young (RY) closure [44]. Computing thus g(r), we obtain S(0) from the integration in equation (11) supplemented with finite size scaling, and the pressure from To assess the validity of both MC and integral equation routes, we show in figure B.1 the results of our S(0) calculations in which the PB-cell effective parameters have been taken in the Yukawa potential. The results of the MC simulations compare favourably with the solutions of the OZ equation for all values of κa considered. Equally good agreement between the results of both methods was found for all other curves presented in this work, but in order not to overload the graphs, we have shown only the OZ results.

Appendix C. Integral equations theory
The structure of a fluid can be obtained by experimental techniques, computer simulations, or by solving numerically the OZ equation. The inhomogeneous OZ equation for a mono-component fluid is given by [43] h(r 1 , r 2 ) = c(r 1 , r 2 ) + where h(r 1 , r 2 ) and c(r 1 , r 2 ) are the total and direct correlation functions, respectively, between a particle located at r 1 and a particle located at r 2 . ρ(r 3 ) is the local density of particles in the system. In an homogeneous and isotropic system the total and direct correlation functions depend only on the relative distance between particles and the local density takes the average value ρ(r 3 ) = ρ, where ρ is the mean density of particles. Then, equation (C.1) reduces to The total correlation function is related to the local structure of the system by means of the relation h(r) = g(r) − 1, where g(r) is the radial distribution function. Also, h(r) is connected with the structure factor, S(q), through the relation S (q) = 1 +h(q), whereh(q) is the Fourier transform of h(r). The structure factor can be measured, for example, by light scattering experiments. Equation (C.2) is an integral equation with two functions, h(r) and c(r), which needs an additional relation between the total and direct correlation functions to close the set of equations. The general closure relation for the equation (C.2) is given by [43]  where βu(r) is the pair potential between particles and the function B(r) is the so-called bridge function, which depends on particle density and, in general, is unknown. There are many approximations to the bridge function, such as Percus-Yevick (PY), hyper-netted chain (HNC) and RY closure relations [43,44]. It is well-known that the HNC and RY relations work well when the interaction between particles is only repulsive beyond the hard-core interaction. The RY approximation [44] is given by where the function f(r) = 1 − exp(−αr). The RY closure is a mixture between the PY and HNC closure relations. For example, when α = 0 equation (C.5) reduces to PY approximation. As α increases, f(r) approaches 1, and equation (C.4) reduces to HNC approximation. The mixture parameter, α, is fixed by demanding that the isothermal compressibility from the virial route, 1 − ρc(q = 0) = S(q = 0), and from the compressibility route (see equation (3)), are both equal.
In the case of density-dependent potential of interaction, as is the case here, it is important to compute the density derivatives leading to the above compressibilities at constant potential of interaction.
For a given density ρ the integral equation (C.2), together with the RY closure relation, is numerically solved by converting the OZ equation in an algebraic equation. This is done by taking the Fourier transform of equation (C.2), h(q) =c (q) 1 − ρc(q) . (C.5) Then, by applying the inverse Fourier transform to equation (C.5) we get the desired solution. However, for strongly interacting systems, as in our case, the direct application of equation (C.5) gives noisy solutions. Then, we divide the density ρ in small steps of size ρ. For each given sub-density we solve iteratively the OZ by using a five-parameter version of the Ng-method [45] until the desired density is reached. At each step of the iteration, the pair distribution function is determined and used together with the RY closure relation in order to compute the new direct correlation function, c(r). To ensure rapid convergence, the value of c(r) at the previous step is taken as an initial guess.