Screening of Spherical Colloids beyond Mean Field -- A Local Density Functional Approach

We study the counterion distribution around a spherical macroion and its osmotic pressure in the framework of the recently developed Debye-H"uckel-Hole-Cavity (DHHC) theory. This is a local density functional approach which incorporates correlations into Poisson-Boltzmann theory by adding a free energy correction based on the One Component Plasma. We compare the predictions for ion distribution and osmotic pressure obtained by the full theory and by its zero temperature limit with Monte Carlo simulations. They agree excellently for weakly developed correlations and give the correct trend for stronger ones. In all investigated cases the DHHC theory and its computationally simpler zero temperature limit yield better results than the Poisson-Boltzmann theory.


I. INTRODUCTION
The screening of charged macromolecules in an electrolyte solution is a long standing problem which has prompted many attempts aiming at a theoretical explanation. In their pioneering work Gouy [1] and Chapman [2] used what is now referred to as Poisson-Boltzmann (PB) theory as the basis for a mean field treatment of the electrical double layer. This approach found its culmination about thirty years later in the famous DLVO theory of charged colloids [3,4]. The major flaw of these mean field approaches is their neglect of correlations between the ions. The first attempt to work out such correlations for homogeneous electrolytes are due to Debye and Hückel [5], whose work remarkably (and at first glance confusingly) is also based on (linearized) Poisson-Boltzmann theory. In the inhomogeneous case integral equation theories [6,7,8,9] and recently field theories [10] have become very popular in calculating correlation corrections to mean field double layers. However, in order to make progress and calculate physical quantities, approximations have to be made which, unfortunately, instead of clarifying the physics sometimes tend to obscure it. Moreover, since in some of these methods, the free energy is not defined in a unique way, it becomes impossible to determine the specific role played by each source of correlations in the system. It would therefore be desirable to have a theoretical framework which retains the simplicity of the early attempts, but also accommodates correlation effects. This is the case for density functional theories. It is possible to rigorously rewrite the partition function of, say, a system of charged colloids, as a density functional [11], in which the contribution beyond mean field is seen to be expressible as an additive correlation correction to the free energy density, whose functional form is of course unknown and for which one has to make a reasonable ansatz. The spirit is very similar to the fundamental problem of integral equations, where one also has to make an educated guess (namely, the closure relation), but in the functional case the ansatz involves a free energy density rather than a relation between two-and three-point functions. It thus relies on a different kind of intuition and thus permits complementary insight.
One suggestion for such a functional correction has been made by Nordholm [12]. It relies on a Debye-Hückel treatment of the One Component Plasma (OCP) [13,14,15], in which the short-distance failure of linearization is cleverly overcome by postulating a correlation hole. Since beyond a certain density the resulting OCP free energy density is a concave function of density, this favors the development of inhomogeneities. In the pure OCP these are balanced by the homogeneously charged background. However, if one uses the OCP free energy density as a correlation correction to the mean field functional describing the double layer at a charged surface, one has all the charge opposite to the counterions located on that surface, rather than homogeneously distributed as a stabilizing background. The consequence is that the double layer becomes unstable and all ions collapse onto the surface, an effect which has been termed "structuring catastrophe" [16,17].
To circumvent this instability without losing the physical transparency of a local functional, we recently proposed the Debye-Hückel-Hole-Cavity (DHHC) theory [18], in which we suggested a convex correlation functional. This was achieved by excluding the homogeneous background from a region of radius a around the central ion during the Debye charging process. For counterions with size we identified a tentatively as the ion diameter. We then applied our theory to the screening of a charged rod by its counterions. Comparisons of the ionic charge distribution obtained showed a very good agreement with the simulations for both monovalent and trivalent counterions.
In this paper we test our theory for a different geometry: charged spherical colloids with point-like counterions. In general, colloidal systems exhibit a rich phase behavior. The particles can agglomerate at high densities, generally an irreversible process, but they may also show a reversible liquid-vapor phase separation similar to the one present in simple molecular liquids. In order to prevent them from simply falling out of solution, one needs some kind of repulsion between the particles. Introducing charged groups at the surface of the colloid is one way to do that. The large gain in entropy following the dissociation of a vast number of counterions into solution stabilizes the system, because an aggregation of colloids into a small sub-volume would -for reasons of global charge neutrality -also require the counterions to occupy this small volume and thereby give up much entropy. Of course, the final state of the system is always a balance between energy and entropy, and if electrostatic interactions are strong, they will ultimately overcome entropy and lead to aggregation of the colloids [19,20,21,22,23]. The resulting phenomenon of "like charge attraction" has received much attention, but it is of course only mysterious if one forgets that the entire system is neutral. Admittedly, confusion persists about whether such a phase separation could also happen within mean field theory. Even though rigorous proofs exist that PB theory will not permit attraction between like charged macroions under reasonably general circumstances [24,25,26], and that in a cell model treatment the compressibility will be positive [27], it has been claimed that an expansion of the free energy of a charged colloidal suspension into zero-, one-, two-etc. body terms will contain configuration independent volume terms, which may drive a phase separation even though the pair terms are purely repulsive [28,29,30]. Since unfortunately all these derivations rely on a linearization of PB theory, which might render the findings as artifacts [31,32,33,34], the issue appears to be open yet.
All these phenomena ultimately depend on the screening produced by the ionic cloud, which in turn depends on the geometry of the system. In this regard, a charged spherical colloid differs from a charged rod in two fundamental ways: the electrostatic potential and the spatial extension. The logarithmic potential present in the case of charged cylinders leads to the phenomenon known as Manning condensation [35,36]. If the line charge density exceeds a critical threshold, a certain fraction will remain loosely associated with the rod, even at infinite dilution, and renormalize the rod charge. A quantitative PB treatment of this provides a unique criterion for defining the effective charge of the system, even at finite densities [37,38].
The situation is different for charged spherical colloids, which lose all their counterions in the limit of infinite dilution; thus, the colloidal charge does not get renormalized. Still, on often talks about effective charges, which mimic the stronger condensation of nonlinear theory within a linearized treatment [39,40,41,42,43,44]. That, however, is clearly not a physical but rather a formal renormalization, necessitated by the simplified linear treatment, and is thus a different story.
Another important difference between the spherical and the cylindrical symmetry lies in the spatial extend. If a charged rod is infinitely long (as is usually assumed in theoretical treatments), the number of counterions at any given distance from the rod is always infinite. In contrast, for a charged spherical colloid the number of counterions at any distance is always finite, since of course there is no direction along which the colloid is infinite. Hence, fluctuations of the radial charge density are more likely to be important in the spherical case.
The systems we will consider here are strongly charged colloids with point-like ions of some specific valence and no added salt inside a spherical cell. Since all the particles are limited to be within one cell, correlations between different macroions and between microions belonging to different cells are not present. In our treatment we will thus exclusively focus on questions regarding the description of a single double layer. Furthermore, for point-like ions the interpretation of our cutoff parameter, a, can obviously no longer be the particle diameter. We will introduce an alternative prescription for a, based again on local density considerations and keeping in mind that its entire purpose is to prevent the functional from becoming unstable.
We also derive an approximated version of our correlation functional, namely, its zero temperature limit. It has the huge advantage that it can be calculated analytically, while still predicting ion profiles quite close to the full DHHC expression for a wide range of parameters. It also demonstrates the spirit of our stabilization correction very directly.
Finally, we compare our predictions for ion profiles with Monte Carlo simulations, in which we independently vary valence v and plasma parameter Γ 2d = πσℓ 2 B v 3 , where σ is the density of surface charges and ℓ B is the Bjerrum length. It has been shown that beyond Γ 2d ≃ 2.26 the force-distance curves between charged plates cease to be monotonic, and beyond Γ 2d ≃ 2.45 attractions set in [45]. These effects result from correlations between different double layers (like, for instance, ion interlocking [46,47]), which we cannot account for, and it has in fact been shown that they cannot be described within a local density functional theory with a convex correlation correction [26]. However, for the description of a single double layer the regime of applicability of our theory is larger, even though it clearly must fail for too high coupling.
The paper is organized as follows. In Sec. II the DHHC correlation functional is revisited and its zero temperature limit is introduced. It is then applied as a local correlation correction to the problem of screening of charged colloids in Sec. III. The case of point-like ions is discussed in detail and the new expression for a is proposed. Technical details of the simulations are described in Sec. IV. The results of the simulations, full theory and zero temperature limit are compared in Sec. V, and we end with our conclusions Sec. VI.

II. THE DEBYE-HÜCKEL HOLE-CAVITY (DHHC) THEORY REVISITED
The one component plasma consists of N identical point-particles of valence v and (positive) unit charge q inside a volume V with a uniform neutralizing background of charge density −vqn B and dielectric constant ε. As a first approximation the free energy of this system can be derived in the framework of the Debye-Hückel approach. Then, the electrostatic potential ψ created by some ion, fixed at the origin for instance, and all its surrounding ions satisfies the spherically symmetric Poisson equation ∇ 2 ψ(r) = ψ ′′ (r) + 2 r ψ ′ (r) = −4πρ(r)/ε. The charge density has a contribution from the central ion, vqδ(r), a contribution from the surrounding ions which are distributed -within mean field theory! -according to the Boltzmann factor n PB (r) = vqn B exp{−βvqψ(r)}, and finally from the charged background. Inserting this into the Poisson equation and linearizing the exponential yields the linearized Poisson-Boltzmann equation The solution of Eqn. (1) is the well known expression ψ(r) = vq e −κr /εr. However, the problem with Debye-Hückel theory is that the condition for linearization is obviously not satisfied for small r, where the potential is large. Indeed, since all ions have the same sign of charge, this implies that the particle density becomes negative and finally diverges at the origin. This defect was overcome by the Debye-Hückel-Hole theory [12], which artificially postulates a correlation hole of radius h around the central ion into which no other ions are allowed to penetrate. In this case the charge density is given by The hole size h is fixed by excluding particles from a region where their Coulomb energy is larger than k B T , which gives 1 + κh = (1 + 3κℓ) 1/3 . Once the potential at the position of the central ion is known, the electrostatic contribution to the free energy density, f DHH (n), can be obtained by the Debye charging process [5,17].
The simple Debye-Hückel-Hole analysis of the onecomponent plasma theory offers considerable insight into ionic systems and is in good agreement with Monte-Carlo simulations [48] when fluctuations on the charge density are not relevant [49]. In principle, one can attempt to include such fluctuations by applying the bulk densityfunctional theory in a local way. The basic idea is to obtain the density distribution via functional minimization of the free energy The first part, the PB free energy contains the entropy of the mobile ions, the interaction of the small ions with the macroion potential and the mean-field interaction between the counterions. Here V p is the particle volume. The expression f corr in Eqn. (3) accounts for the correlation between the mobile ions. The ion distribution can be derived by minimization of Eqn. (3) under the constraint of charge neutrality. Unfortunately, this variational process does not lead to a well defined density profile if one uses f DHH (n) as the correlation correction f corr . The reason is that f DHH (n) is a concave function beyond n ⋆ ≈ 7.86/ℓ 3 and asymptotically behaves as −n 4/3 . Since there is no stabilizing homogeneously charged background but rather a concentration of opposite charge on the macroion surface, this favors the development of a distribution in which all the ions sit on the surface of the macroion.
The instabilities present in the DHH approach can be properly overcome by recognizing that the failure of this model is due to the too strong requirement of local charge neutrality imposed by the local density approximation: A local fluctuation leading to an increase of particle density implies a corresponding increase in background density. Therefore, the fluctuation is not suppressed by an increase in repulsive Coulomb interactions but quite on the contrary favored by its decrease. To circumvent the instabilities occurring at high densities, we proposed recently a simple solution in which one excludes the neutralizing background from a cavity of radius a placed around the central ion (for details of the derivation of the model see Ref. [18]). In this case, the charge density can be split in three different regions, namely where the hole size h is chosen such as to yield the same screening (i.e., the same amount of charge within h) as the DHH theory, which results in with ω = (1 + 3κℓ) 1/3 . Using this prescription for h, the free energy is obtained by Debye-charging the fluid: where Since the DHHC free energy is a convex function of density, f DHHC can thus be used to account for correlations within a local density approximation.

A. The Zero Temperature Limit DHHC (0)
The fact that the integral in Eqn. (7) has to be solved numerically obstructs a direct view on how thermodynamic stability is actually restored. Luckily, the crucial point can already be seen by focusing on the limit of zero temperature. In this case Eqn. (6) gives the expression for the correlation hole of the DHHC theory. This conveniently implies the potential to vanish outside h. In other words, the region a < r < h contains the right amount of background charge to exactly neutralize the central ion, and it is appropriate to refer to this limit as "complete screening". The potential in the two other regions then simplifies considerably: with the dimensionless scaled density n B given by n B = 4 3 πa 3 n B . After the Debye charging process one obtains the following closed expression for the excess free energy density: Note that the limits a → 0 and n B → ∞ do not commute: For high densities, βf and this concave scaling with density prevents it from being used within a local density approximation. The zero temperature limit thus demonstrates in a clear way the key role played by the cavity of size a, which excludes the uniform background from the vicinity of the central ion.

B. How to choose a proper value for a
Before applying this strategy to various valences and ionic strengths, we need to specify the parameter a. If the counterions have a diameter d, no other charge should be found at a distance r < d. Therefore, in the spirit of the Debye-Hückel theory, we tentatively interpreted a in Ref. [18] as the ion diameter. This choice has led to an excellent agreement with simulations when applied to rod-like polyelectrolytes [18] with mono-, di-, and trivalent counterions (and no added salt), but it is of course infeasible for point ions. In the following we suggest an alternative way to choose a value for a which is independent of excluded volume arguments, and show that this choice yields a good description of our Monte Carlo results and trends.
We already mentioned the crucial role played by a in maintaining the free energy convex. We also have seen in our discussion of the zero-temperature case that this is achieved because a in Eqn. (9) balances the length (3/4πn B ) 1/3 , which is basically the mean distance between ions. One could thus try to selfconsis-tently choose a proportional to the local ion distance, but this would be unsuccessful: The balance would not work, since each density increase would shrink a proportionally, and the collapse could not be stopped. One thus needs a length which is local and somehow related to the ion density -but which does not change as the local ion density changes. This suggests to pick the average distance between ions as predicted by PB theory: a = (3/4πn PB (r)) 1/3 . Our density functional then quite naturally emerges as a next order correction to the mean field result.
After these general considerations on the cutoff a, let us continue with a practical remark. Far away from the charged surface the ion density is always quite low, correlations are weakly developed, and the precise value of a is immaterial. In fact, we only ever need a stabilizing cutoff close to the charged surface, where the ion density is largest. This suggests the following simplification: Instead of using a cutoff function a(n PB (r)) depending on the local PB density, we pick a constant a from a worstcase scenario, namely, the value which it has at contact. This then finally yields the following prescription for a: In fact, since the cutoff will become important in the regime of strong correlations, we could even replace the contact density n PB (r 0 ) by its limiting value 2πℓ B σ 2 , where σ is the density of surface charges [50]. We then find a ℓ strong coupling where is the 2d plasma coupling parameter [15]. Formula (14) nicely demonstrates that in this limit the cavity size, measured in the appropriate length scale ℓ (see also the scaling discussion in the Appendix), is simply another measure of the coupling strength.

III. APPLICATION TO THE SPHERICAL CELL MODEL
Charged spherical colloids are common and well characterizable systems for studying many electrostatic phenomena in pure culture, and they can often serve as simplified models for more complicated systems like polyelectrolytes or proteins. Solutions containing such charged structures are indeed complicated to describe due to the long-range nature of the Coulombic interactions. However, as long as these long range forces are repulsive, the colloids will create large correlation holes ("cells") around themselves which are void of other colloids. In a first approximation one can then decouple the macroion interactions and concentrate on what's going on within a single correlation hole-an approach which is termed "cell model" [51]. The cell picture is known to give a good approximation for many realistic systems, and most of the physics of the system is determined by the screening of the macroion by the microions inside a cell. As a test case for our theory we shall therefore consider a charged spherical colloid of radius r 0 containing Z charged groups, which are neutralized by point-like ions of valence v. This macroion is embedded in the center of a spherical cell of radius R, corresponding thus to a volume fraction φ = (r 0 /R) 3 of colloids.
The thermodynamic behavior of the colloidal system is determined by the distribution of mobile ions around the macroion. This distribution is obtained by minimization of the free energy functional, Eqn. (3). For the colloidal system, the interaction of the small ions with the macroion and the mean-field interaction between the counterions are given by where ψ(r) is the total electrostatic potential at position r and ψ fix (r) = −Zq/εr is the potential due to the charged macroion alone. The inter-particle correlations are taken into account by employing f corr = f DHHC . The minimization itself is accomplished by numerically solving the corresponding Euler-Lagrange-equation. Special care had to be taken to obtain a sufficient accuracy of the rapidly varying density profiles close to the colloid surface. In the following we will concentrate on two observables. The first is the integrated fraction of ions within a radial distance r from the colloid center, which is given by We will sometimes plot P as a function of −1/r, which is the Green function of the spherical Laplacian. This will visually expand the region close to the colloid, but it also has practical advantages when estimating the amount of closely associated ions, see e.g. Ref. [40,52]. Measuring all lengths in the full partition function of the cell model in units of ℓ = ℓ B v 2 reveals that the distribution function P (r) is invariant under a rescaling which keeps the number of counterions N = Z/v, the reduced colloid size r 0 /ℓ, and the volume fraction φ = (r 0 /R) 3 constant (see Appendix). The same holds for PB theory, and it is also true for DHHC theory. In the latter case this not only relies on the form of the DHHC free energy correction (7), but also on our particular choice of a. This invariance property is thus a further support for Eqn. (13).
The second observable we look at is the osmotic pressure Π. For PB like free energy functionals with an additional density term -like our f corr -it is given by [27] For the PB case, f corr ≡ 0, this reduces to the well known fact that the pressure is given by the boundary density [53]. Since this result actually holds rigorously for the full restricted primitive model [50], one could also argue that DHHC theory is an approximate way to calculate the boundary density, and then calculate the pressure from βΠ = n(R), i.e., leave out the additional term nf ′ − f . This would lead to a different result, reminding us that selfconsistency and consistency with other rigorous results cannot generally be achieved. We will always use the internally consistent equation (18) for our pressure calculations.

IV. SIMULATIONAL DETAILS
The systems we study consist of a spherical macroion of radius r 0 and (negative) central charge −Zq. Electroneutrality is ensured by the presence of N = Z/v point-like counterions of valence v, confined inside an impermeable spherical cell of radius R. This also fixes the colloid volume fraction to φ = (r 0 /R) 3 . No additional salt is added. The dielectric constant ε is assumed to be uniform throughout the system, such that no image forces [54] occur. Our choices for the system parameters can be found in Table I Standard canonical MC simulations following the Metropolis scheme [55] were employed to sample the ion distributions. After an initial equilibration time of 200 000 MC steps, where we attempted to move every ion once to a new position, we sampled the system for  The parameters of the simulated systems. The volume fraction was always chosen as φ = (r0/R) 3 = 0.8%. The solid curve is the result of the MC simulation, while the dash-dotted curve is the prediction from PB theory. The inset shows the local density n(r). The increase in the counterion condensation due to correlations is well captured by the DHHC theory (dashed curve) and its zero temperature limit DHHC (0) (dotted curve). The difference in n(r) between the latter two is invisible on the chosen scale, and only DHHC is shown.
1.3 − 2 × 10 6 MC steps, producing 1300-2000 configurations for analysis. We will measure energies in units of k B T and use the coupling length ℓ = ℓ B v 2 as our unit of length (for monovalent ions under aqueous conditions and room temperature we would have ℓ = 7.14Å, and the unit of concentration becomes ℓ −3 = 4.56 M). In the following we will present MC results for the integrated ion distribution, Eqn. (17), and for the pressure, Eqns. (19) and (20).

V. COMPARISON BETWEEN SIMULATIONS AND DHHC THEORY
A. Ion distribution functions Figure 1 shows the integrated charge distribution P (r) for system 1 from Tab. I. The solid curve is the result from the MC simulation, and it lies distinctly above the PB result (dash-dotted curve), indicating a stronger con- Counterion distribution function P (r) for systems 2, 3, and 4 from Tab. I. The line styles are the same as in Figure 1, the counterions are monovalent, and the value of the plasma-parameter Γ 2d is indicated. densation of ions due to correlations neglected in PB theory. Most of this enhancement of ion localization close to the colloid is captured by DHHC theory (dashed curve) or its zero temperature limit DHHC (0) . This is also evident from the local density n(r), which relatively to PB is enhanced at close proximity to the colloid, while it drops below PB at the outer cell boundary. From what we have said in section III this also indicates that the pressure will be lower, and this is indeed what we shall find (see below).
It should be noted that the 2d plasma parameter Γ 2d = 2.5 is already slightly beyond the point where attractions between two planes would arise [45]. We should not expect DHHC theory to work for significantly higher plasma parameters, since it cannot account for effects like attractions [26]. However, we want to point out that here we only aim at properties of a single electrostatic double layer and not at phenomena arising from the interaction between two of them, and in fact the agreement seen in Fig. 1 is very encouraging. It is also quite pleasing that the significantly simpler zero temperature limit DHHC (0) from Eqn. (11) yields essentially the same result as the full DHHC theory.
Due to the scaling invariance of the partition function discussed in the Appendix, a system with e.g. divalent ions and Z = 200 or trivalent ions and Z = 300 (and properly rescaled Bjerrum lengths ℓ B → ℓ B /v 2 ) shows exactly the same distribution function (not shown).
In Figure 2 we show distribution functions P (r) for systems 2-4 from Tab. I. These have monovalent counterions and only differ in their value of the plasma parameter Γ 2d . Clearly, a larger plasma-parameter leads to an increased condensation (the curves are shifted up)-an effect which naturally is already present in PB theory. However, apart from this, at a larger plasma parameter the influence of correlations becomes more important,  Figure 1, the plasma parameter is Γ 2d = 2, and the value of the counterion valence is indicated. For clarity, the PB curve is only shown for v = 3. and therefore the deviation between the PB prediction and the MC result increases for increasing Γ 2d , which is also clearly seen in Fig. 2. Again, this effect is well captured by DHHC theory, which is always much closer to the MC data than to the PB result, even though its accuracy diminishes as Γ 2d becomes large.
In Fig. 3 we show a "complementary" scan, in which we fixed the value of the plasma-parameter Γ 2d = 2, but changed the counterion valence (systems 3, 5, and 6 from Tab. I). Maybe surprisingly, an increase in valence leads to a decrease in condensation if it happens at constant plasma-parameter and colloid charge. If we had changed v from 1 to 2 and simultaneously replaced ℓ B → ℓ B /4 and Z → 2Z, the plasma parameter would also have remained unchanged, but due to the scaling property of the partition function that would actually have been true for the whole distribution function. Instead, we have reduced ℓ B → ℓ B /2 3/2 ≈ ℓ B /2.83 (i.e., a little less strongly), but have failed to increase Z. The net result is that condensation drops slightly. However, since the plasma parameter, which is the best indicator for the strength of correlations, remains the same, the deviation between PB theory and MC simulation are always about the same (not shown in the Figure). And as a consequence, the deviation between DHHC theory, which approximately accounts for correlations, and the MC simulation, which captures them all, is about the same in all three cases, and actually not very big.

B. Osmotic Pressure
Another strategy to check how successful our approach captures correlations is to compute the osmotic pressure. In real systems this pressure will depend on correlations between ions of different cells, something which neither our theory nor actually our simulation (of a single colloid!) takes into account. So in the following by "pressure" we do not, strictly speaking, refer to the bulk pressure of a colloidal suspension at some given volume fraction, but only to the pressure exerted on the rigid wall at r = R of our cell model.
Within the simulations, the pressure is given by the contact density at r = R, which was obtained by fitting the MC density profile close to the cell boundary to a quadratic expression n(r) = c 1 + c 2 (r − R) 2 . An example for how the simulated densities compare to the PB approximation, our analytic DHHC approach, and its simpler zero temperature limit DHHC (0) , can be found in Fig. 4. Table II shows the predictions for the pressure in the case of point-like ions given by PB-, DHHC-, and DHHC (0) -theory, as well as by MC simulations. In the case of PB-theory and MC the pressure is simply the density at the outer boundary, while for the DHHC approach we employ Eqn. (19) and for DHHC (0) Eqn. (21). These data, as well as Fig. 4, demonstrate that -as anticipated -the simulated pressures lie below the PB prediction. This decrease in pressure is rather accurately captured by our functional f DHHC and by its zero temperature limit, f (0) DHHC , with the MC result lying significantly below PB and (in these cases) below DHHC and above DHHC (0) . The difference between the two correlation corrected approaches is consistent with the idea that entropic effects neglected in f (0) DHHC would push ions away from the macroions, or in other words, that the zero temperature limit implies stronger correlations than the regular DHHC theory and therefore yields even lower pressures.

VI. CONCLUSIONS
In this paper we showed how to apply our previously proposed local density functional approach based on a stable correlation correction to a spherical macroion confined in a spherical cell. One of the crucial parameters in this theory is the size a of the exclusion cavity of the background charge density. For point-like ions, we suggest to associate the exclusion region with the mean distance between ions as predicted by PB theory, and for simplicity use the value present at colloidal contact.
By going to the zero temperature limit we were able to derive an even simpler free energy functional F (0) DHHC , which is almost as good as the full DHHC theory, but much easier to handle. We also derived exact expressions for the osmotic pressure in this system. We successfully compared our predictions to simulations of the same model and compared the integrated counterion density and the osmotic pressure values for two complementary "scans" of the coupling strength, namely valence and plasma parameter. We demonstrated that our local density functional approach based on a stable correlation correction leads to a major improvement over the PB prediction.
The canonical partition function Z of the colloid surrounded by its counterion is given by: where N = Z/v is the total number of counterions and the Hamiltonian H = T +V splits into kinetic and potential degrees of freedom. In the classical description employed here the kinetic part T will contribute the usual factor λ −3N to the partition function, where λ is the thermal de Broglie wavelength. The potential energy can be expressed as  After rescaling all length by ℓ, i.e. introducing x := r/ℓ, the total partition function can be written as In this form it becomes evident that appropriately scaled thermal observables like the integrated charge density (measured in units of ℓ −3 ) or the pressure (measured in units of k B T ℓ −3 ) are invariant under system changes which leave the number of counterions N , the rescaled colloid size x 0 = r 0 /ℓ, and the volume fraction φ fixed.
Poisson-Boltzmann theory shows the same invariance property, as does the approximate density functional theory we are proposing in this paper.