# Density-functional theory for attraction between like-charged plates

## I Introduction

Solutions containing macromolecules are ubiquitous in the everyday life. From food colloids to the DNA, we are surrounded by these giant molecules which directly or indirectly govern every aspect of our lives. In many cases the macromolecules in solution posses a net charge. The electrostatic repulsion between the polyions is, often, essential to stabilization of colloidal suspensions. In the biological realm the electrostatics is responsible for the condensation of the DNA and formation of actin bundles, while various physiological mechanisms depend on the electrostatic interactions between the proteins and the microions. In spite of their ubiquity, our understanding of polyelectrolyte solutions is far from complete.

The effort to fathom the role of electrostatics as it applies to the colloidal suspensions goes back over half a century to the classic works of Derjaguin and Landau  and of Verwey and Overbeek (DLVO). These in turn were based on the pioneering studies of Gouy  and Chapman  of double layers in metal electrodes. Following these early contributions, a large effort has been devoted to solve the Poisson-Boltzmann (PB) equation in various geometries. The mean-field treatment, based on the solution of the PB equation, suggests that the interaction between two equally charged macroions in a suspension containing counterions is always repulsive. In recent years, however, this dogma began to be questioned based on simulations, analytical calculations  and experiments, which indicated that for small distances and large charge densities, two like-charged polyions might actually attract!

The fundamental goal of this paper is to demonstrate that this attraction is linked to the correlations between the microions omitted in the mean-field theories, and to establish the conditions under which the attraction becomes possible. We shall consider the interaction between two infinite uniformly charged plates confining their own point-like counterions. The mean-field approximation for this system is obtained by solving the PB equation which, due to the planar symmetry, can be done analytically. Once the density profile is obtained, all the other thermodynamic quantities can be easily derived. Thus, it is not difficult to demonstrate that the pressure at the mean-field level, in units of energy, is simply the density of counterions at the mid-plane between the plates. Since this is always positive, no attraction is possible within the mean-field theory.

Realization that the correlations between the counterions can strongly modify the mean-field predictions goes back a number of years. One of the first approaches proposed by Kjellander and Marčelja  was to include the correlations through the numerical solution of the Anisotropic Hypernetted Chain Equation (AHNC). These authors found that the force per unit area (pressure) can become negative in the presence of divalent counterions. Monte Carlo (MC) simulations performed by Guldbrand et al.  also indicate that as the surface charge density is increased, the pressure decreases if the distance between the charged surfaces is sufficiently small. As in the case of the AHNC calculations, attraction was found only in the presence of divalent counterions. These authors, however, did not analyze the case of very high charge density and short distance between the plates. In addition, since in the above calculations it is difficult to separate the different physical contributions to the pressure, the mechanism that drives the attraction remains unclear.

A different theoretical approach which attempted to shed some light on the mechanism of attraction was advanced by Stevens and Robbins. These authors proposed a density-functional theory similar to the one often employed in studies of simple liquids. This approach introduces a grand-potential free energy, $`\mathrm{\Omega }\left[\rho (𝒓)\right]`$, which is a functional of the non-uniform density of counterions $`\rho (𝒓)`$. The equilibrium properties of the system are obtained through the minimization of the total free energy. The practical problem with this method is that the exact form of the functional is not known. When the correlations between the microions are omitted, the minimization of the grand potential, $`\mathrm{\Omega }_{\mathrm{PB}}`$, becomes trivial and leads to the usual PB equation. In order to account for the correlations between the counterions, Stevens and Robbins  appealed to the Local Density Approximation (LDA). Within this approach an additional contribution, $`f_{\mathrm{LDA}}`$, is added to the mean-field expression, $`\mathrm{\Omega }_{\mathrm{PB}}`$. The expression for $`f_{\mathrm{LDA}}`$ adopted by Stevens and Robbins was obtained through the extrapolation of the MC data for the homogeneous One-Component Plasma (OCP), but with the homogeneous density replaced by an inhomogeneous density profile. The minimization of the free-energy functional allowed them to determine the density profile, $`\rho (𝒓)`$, and the pressure, $`P_{\mathrm{OCP}}`$. The LDA, however, is not without its own problems. The major drawback of this approach is that, for short distances and high charge densities, the LDA is unstable. The reason for the instability is due to the fact that as the density of counterions in the vicinity of the plates increases, the chemical potential decreases, what attracts more particles to the region. This, in turn, leads to an unphysical “chain reaction” where all the counterions condense onto the plates. Clearly, when the distance between the counterions becomes smaller than some threshold value, $`s_{\mathrm{corr}}`$, the LDA ceases to be a reliable approximation.

An improvement over the LDA is, the so called, Weighted Density Approximation (WDA). In this case, the excess free energy is taken to be a function of an average density, $`\rho _w(𝒓)=\mathrm{d}^3𝒓^{}w(|𝒓𝒓^{}|)\rho (𝒓^{})`$, averaged over a region of radius $`s=s_{\mathrm{corr}}`$, where the interactions between the counterions are the strongest. The difficulty in the practical implementation of this scheme is the determination of a proper weight function. The simplest possible form for $`w(|𝒓𝒓^{}|)`$, used by Stevens and Robbins, was to assume that this function has a long-range variation comparable to the wall separation. In this case, the weighted density $`\rho _w(𝒓)`$ is approximated by the homogeneous density independent of $`𝒓`$. However, when the walls are not close, $`L>s_{\mathrm{corr}}`$, the weighted function is no longer uniform and the approximation adopted by Stevens and Robbins becomes unrealistic.

A beautiful explanation of the attraction between like-charged plates has been recently advanced by Rouzina and Bloomfield. These authors present a picture of attraction as arising from the ground-state configuration of the counterions. Clearly at zero temperature the counterions will recondense onto the surface of the plates forming two intercalating Wigner crystals. The authors advance a hypothesis that even at finite temperatures, relevant to the common experimental conditions, the attraction is still governed by the zero-temperature correlations. A somewhat different formulations based on field-theoretic methodology have also been proposed. In these approximations the attraction arises as a result of correlated fluctuations in the counterion charge densities. Although providing a nice qualitative explanation of the origin of the attraction, these simple theories fail to yield a quantitative agreement with the simulations.

In this paper we propose a different form of the weighted-density approach, which rectifies the problems of the earlier theories while still remaining numerically tractable. The excess free energy and the weight function, $`w(|𝒓𝒓^{}|)`$, are derived from the Debye-Hückel-Hole (DHH) theory of the OCP. The density profile is determined by minimizing the free-energy density with respect to the local density. Once the density profile is obtained, the free energy of the system is calculated by inserting it into the expression for the free-energy functional. Given the free energy, all the thermodynamic properties of the system can be easily calculated. A careful analysis of the behavior of the pressure as a function of the charge density and the distance between the plates allows us to explore the nature and the origin of the attraction.

The remainder of the paper is organized as follows. The model and the PB approximation for the density-functional approach are described in Sec. II. The WDA is introduced and applied in Sec. III. Our results and conclusions are summarized in Sec. IV.

## II The Poisson-Boltzmann Approach

We consider two large, charged, thin surfaces each of area $`𝒜`$, separated by a distance $`L`$ (see Fig. 1). The two plates with a negative surface charge density, $`\sigma `$, confine positive point-like monovalent counterions with charge $`e`$. The overall charge neutrality of the system is guaranteed by the constraint

$$_{L/2}^{L/2}dz\rho (z)=\frac{2\sigma }{e},$$
(1)   
where $`\rho (z)`$ is the local number density of counterions and $`z`$ is the Cartesian coordinate perpendicular to the plates. The space between the plates is assumed to be a dielectric continuum of constant $`\epsilon `$.

In order to explore the thermodynamic properties of the system, we use a density-functional approach. The grand potential of the system is

$$\mathrm{\Omega }\left[\rho \right][\rho ]\mu N,$$
(2)   
where $`N`$ is the total number of counterions, $`\mu `$ is their chemical potential and the functional $``$ is derived from the free-energy density of the homogeneous system, with the uniform density of counterions, $`\rho _\mathrm{c}=N/L𝒜`$, replaced by the local density $`\rho (z)`$. For dilute systems, the ionic correlations can be neglected and the grand-potential functional (per unit area) becomes

$``$

$$\frac{\beta \mathrm{\Omega }\left[\rho \right]}{𝒜}=_{L/2}^{L/2}dz\rho (z)\left\{\mathrm{ln}\left[\mathrm{\Lambda }^3\rho (z)\right]1\right\}+\frac{\beta }{2}_{L/2}^{L/2}dz\varphi (z)\left[e\rho (z)+q(z)\right]\beta \mu _{L/2}^{L/2}dz\rho (z),$$
(3)   
$``$

where the electrostatic potential,

$$\varphi (𝒓)=\mathrm{d}^3𝒓^{}\frac{e\rho (𝒓^{})+q(𝒓^{})}{\epsilon |𝒓𝒓^{}|},$$
(4)   
due to the symmetry of the problem, depends only on the $`z`$ coordinate. $`\mathrm{\Lambda }`$ is the de Broglie thermal wavelength of the counterions, $`\beta =1/k_BT`$ and $`q(z)=\sigma \left[\delta (zL/2)+\delta (z+L/2)\right]`$ is the surface charge density of the plates. The functional minimization of this expression,

$$\frac{1}{𝒜}\frac{\delta \beta \mathrm{\Omega }}{\delta \rho (z)}=0,$$
(5)   
produces the optimum density profile,

$$\rho (z)=\rho _0\mathrm{exp}\left[\beta e\varphi (z)\right].$$
(6)   
The constant $`\rho _0`$ is determined from the overall charge-neutrality condition, Eq. (1),

$$\rho _0\frac{2\sigma }{e{\displaystyle _{L/2}^{L/2}}dz\mathrm{exp}\left[\beta e\varphi (z)\right]}.$$
(7)   
The electrostatic potential is obtained by solving the Poisson equation,

$$\frac{\mathrm{d}^2\varphi (z)}{\mathrm{d}z^2}=\frac{4\pi }{\epsilon }\left[e\rho (z)+q(z)\right],$$
(8)   
with the distribution of free ions given by Eq. (6). We find

$$\varphi (z)=\frac{1}{\beta e}\mathrm{ln}\left[\mathrm{cos}^2\left(\frac{zz_0}{\lambda }\right)\right]\varphi _0,$$
(9)   
where $`\varphi _0`$ is the reference potential, which we will set to zero. Here $`\lambda =1/\sqrt{2\pi \lambda _B\rho _0}`$ and $`\lambda _B=\beta e^2/\epsilon `$ is the Bjerrum length. Eq. (8) has to obey two boundary conditions, namely,

$`E(z=0)`$ $`=`$ $`0,`$ (10)   
$`E\left(z=\pm {\displaystyle \frac{L}{2}}\right)`$ $`=`$ $`\pm {\displaystyle \frac{4\pi \sigma }{\epsilon }}.`$ (11)   
From the first equation, the electric field vanishes at the mid-plane and, therefore, $`z_0=0`$. The second equation imposes the discontinuity of the electric field at both charged surfaces, leading to

$$\frac{1}{\lambda }\mathrm{tan}\left(\frac{L}{2\lambda }\right)=\frac{2\pi \sigma \lambda _B}{e}.$$
(12)   
The potential at a point $`z`$ is, then, given by

$$\varphi (z)=\frac{1}{\beta e}\mathrm{ln}\left[\mathrm{cos}^2\left(\frac{z}{\lambda }\right)\right],$$
(13)   
with $`\lambda `$ the root of Eq. (12). The optimum density profile derived from this potential,

$$\rho (z)=\frac{\rho _0}{\mathrm{cos}^2(z/\lambda )},$$
(14)   
can now be substituted into the free-energy functional, allowing the calculation of the total free energy. The thermodynamic properties of the system can be determined from a suitable differentiation of the total free energy. For example, the force between the two plates is given by the minus derivative of the free energy with respect to the separation $`L`$ between the two surfaces. This differentiation leads to a particularly simple expression for the force per unit of area (or pressure),

$$\beta P=\rho _0.$$
(15)   
We note that although it might be tempting to attribute this simple result to the contact theorem, this is not the case, since the conditions under which this theorem holds are violated in the present geometry; Eq. (15) is purely a mean-field result.

## III The Weighted-Density Approximation

For dense systems, the correlations between the microions become relevant. For instance, if a counterion is present at the position $`𝒓`$, due to the electrostatic repulsion, the probability that another counterion is located in its vicinity is drastically reduced. The correlations in the positions of the counterions reduce the mean-field estimate of the electrostatic free energy. No exact method exists for calculating this excess contribution. The simplest approximation, the LDA, consists of adding to the Eq. (3) a local functional,

$$f_{\mathrm{LDA}}=_{L/2}^{L/2}dz\rho (z)f_{\mathrm{corr}}\left[\rho (z)\right],$$
(16)   
where $`f_{\mathrm{corr}}\left[\rho (z)\right]`$ is the correlational free energy per particle. Within the LDA one normally uses the expression derived for the homogeneous system, in which the uniform density $`\rho _\mathrm{c}=N/L𝒜`$ is replaced by the local density profile $`\rho (𝒓)`$. Unfortunately, as was mentioned above, the LDA is unstable when the one-particle density $`\rho (𝒓)`$ is a rapidly varying function of the position. For example, for high surface charge densities, the minimization of the grand potential has no solution. To circumvent this and related problems intrinsic to the LDA, Tarazona  and Curtin and Ashcroft  proposed a WDA, in which the free-energy density, $`f_{\mathrm{LDA}}`$, is replaced by

$$f_{\mathrm{WDA}}=_{L/2}^{L/2}dz\rho (z)f_{\mathrm{corr}}\left[\rho _w(z)\right].$$
(17)   
The fundamental difference between the LDA and the WDA is that the latter is assumed to depend not on the local density $`\rho (𝒓)`$, but on some average density within the neighborhood of the point $`𝒓`$,

$$\rho _w(𝒓)=\mathrm{d}^3𝒓^{}w[|𝒓𝒓^{}|;\rho (𝒓)]\rho (𝒓^{}).$$
(18)   
This provides a control mechanism which prevents an unphysical, singular, buildup of concentration at one point. The grand potential is obtained by adding the excess free energy per area, given by Eq. (17), to the Eq. (3), 

$`{\displaystyle \frac{\beta \mathrm{\Omega }\left[\rho \right]}{𝒜}}`$ $`=`$ $`{\displaystyle _{L/2}^{L/2}}dz\rho (z)\left\{\mathrm{ln}\left[\mathrm{\Lambda }^3\rho (z)\right]1\right\}`$ (22)   
$`+{\displaystyle \frac{\beta }{2}}{\displaystyle _{L/2}^{L/2}}dz\varphi (z)\left[e\rho (z)+q(z)\right]`$   
$`+\beta {\displaystyle _{L/2}^{L/2}}dz\rho (z)f_{\mathrm{corr}}\left[\rho _w(z)\right]`$   
$`\beta \mu {\displaystyle _{L/2}^{L/2}}dz\rho (z).`$   
Minimization of this expression leads to the optimum particle number density, 

$$\rho (z)=\rho _0\mathrm{exp}\left[\beta e\varphi (z)\beta \mu _{\mathrm{ex}}(z)\right],$$
(23)   
where the excess chemical potential derived from $`f_{\mathrm{WDA}}`$, Eq. (17), is

$`\mu _{\mathrm{ex}}(z)`$ $`=`$ $`{\displaystyle \frac{\delta f_{\mathrm{WDA}}}{\delta \rho (z)}}`$ (24)   
$`=`$ $`f_{\mathrm{corr}}\left[\rho _w(z)\right]+{\displaystyle _{L/2}^{L/2}}dz^{}\rho (z^{}){\displaystyle \frac{\delta f_{\mathrm{corr}}\left[\rho _w(z^{})\right]}{\delta \rho (z)}},`$ (25)   
and the normalization coefficient is

$$\rho _0\frac{2\sigma }{e{\displaystyle _{L/2}^{L/2}}dz\mathrm{exp}\left[\beta e\varphi (z)\beta \mu _{\mathrm{ex}}(z)\right]}.$$
(27)   
The electrostatic potential satisfies the Poisson equation, Eq. (8), with the charge density given by the Eq. (23). Integrating the Poisson equation over a rectangular shell of area $`𝒜`$ and width $`z`$, and appealing to the Gauss’ theorem, an integro-differential equation for the electric field $`E(z)`$ can be obtained,

$$(\overline{z})=4\pi \overline{\sigma }\frac{{\displaystyle _0^{\overline{z}}}d\overline{z}^{}\mathrm{exp}\left[\overline{\mu }_{\mathrm{ex}}(\overline{z}^{})+{\displaystyle _0^{\overline{z}^{}}}d\overline{z}^{\prime \prime }(\overline{z}^{\prime \prime })\right]}{{\displaystyle _0^{\overline{L}/2}}d\overline{z}^{}\mathrm{exp}\left[\overline{\mu }_{\mathrm{ex}}(\overline{z}^{})+{\displaystyle _0^{\overline{z}^{}}}d\overline{z}^{\prime \prime }(\overline{z}^{\prime \prime })\right]},$$
(28)   
where $`e\beta \lambda _BE`$, $`\overline{\sigma }\sigma \lambda _B^2/e`$, $`\overline{z}z/\lambda _B`$, $`\overline{L}L/\lambda _B`$ and $`\overline{\mu }_{\mathrm{ex}}\beta \mu _{\mathrm{ex}}`$. The local density $`\rho (z)`$, which enters in the calculation of the excess chemical potential, Eq. (24), can be obtained from the derivative of the electric field, since $`𝑬(𝒓)=4\pi e\rho (𝒓)/\epsilon `$. The Eq. (28) explicitly fulfills the two boundary conditions: $`(0)=0`$ and $`(\pm \overline{L}/2)=\pm 4\pi \overline{\sigma }`$.

The solution of this equation depends on the specific form of the excess free-energy density and the weight function $`w\left(|𝒓𝒓^{}|\right)`$. For the homogeneous OCP the electrostatic free energy can be easily obtained using the DHH theory of Nordholm. This is a simple linear theory based on the ideas of Debye and Hückel. The electrostatic potential of the OCP is assumed to satisfy a linearized PB equation. As a correction for the linearization, Nordholm postulated the existence of an excluded-volume region of size $`s_{\mathrm{corr}}`$, from which all other ions are excluded. The size of this region is such that the electrostatic repulsion between two counterions is comparable to the thermal energy. Recent calculations using a generalized Debye-Hückel theory indicate that this exclusion region is responsible for the oscillations observed in the structure factor of the OCP at high couplings. Following Nordholm, we find

$$s_{\mathrm{corr}}=\frac{1}{\kappa _D}\left(1+3\lambda _B\kappa _D\right)^{1/3}\frac{1}{\kappa _D},$$
(29)   
where $`\kappa _D=\sqrt{4\pi \lambda _B\rho _\mathrm{c}}`$ is the inverse of the Debye length. The excess free energy per particle is calculated to be

$``$

$$\beta f_{\mathrm{OCP}}=\frac{1}{4}\left[1+\frac{2\pi }{3\sqrt{3}}+\mathrm{ln}\left(\frac{\omega ^2+\omega +1}{3}\right)\omega ^2\frac{2}{\sqrt{3}}\mathrm{tan}^1\left(\frac{2\omega +1}{\sqrt{3}}\right)\right],$$
(30)   
$``$

where $`\omega =(1+3\lambda _B\kappa _D)^{1/3}`$. The correlational free energy per particle for the WDA, $`f_{\mathrm{corr}}`$, which appears in (24), is obtained by replacing $`\rho _\mathrm{c}`$ by $`\rho _w(z)`$ in the expression (30), that is, $`f_{\mathrm{corr}}\left[\rho _w(z)\right]=f_{\mathrm{OCP}}\left[\rho _\mathrm{c}\rho _w(z)\right]`$.

To obtain the weighted function  we require that the second functional derivative of the free energy $``$ in the limit of homogeneous densities, 

$`{\displaystyle \frac{\delta ^2\beta }{\delta \rho (𝒓)\delta \rho (𝒓^{})}}`$ $`=`$ $`{\displaystyle \frac{\delta ^3(𝒓𝒓^{})}{\rho (𝒓)}}+w(|𝒓𝒓^{}|){\displaystyle \frac{\delta \beta \mu _{\mathrm{ex}}(𝒓)}{\delta \rho (𝒓^{})}}`$ (32)   
$`+{\displaystyle \frac{\lambda _B}{|𝒓𝒓^{}|}},`$   
produces the direct correlation function $`C_2(𝒓)`$ of the homogeneous system,

$$\frac{\delta ^2\beta }{\delta \rho (𝒓)\delta \rho (𝒓^{})}=\frac{\lambda _B}{|𝒓𝒓^{}|}C_2(𝒓𝒓^{}).$$
(33)   
Following Groot, we find that a reasonable approximation for the weight function is

$$w(r)=w(|𝒓|)=\frac{3}{2\pi s_{\mathrm{corr}}^2}\left(\frac{1}{r}\frac{1}{s_{\mathrm{corr}}}\right)\mathrm{\Theta }(s_{\mathrm{corr}}r),$$
(34)   
where $`\mathrm{\Theta }(x)`$ is the Heaviside step function. It is important to remember that the radius of the excluded-volume region, $`s_{\mathrm{corr}}`$, is now a function of the position, since the average density $`\rho _\mathrm{c}`$, which appears in Eq. (29), is replaced by $`\rho (z)`$, the local density of counterions, see Eq. (18). Taking advantage of the planar symmetry of the system, the expression for the weighted density can be written explicitly as a one-dimensional quadrature,

$``$

$`\rho _w(z)`$ $`=`$ $`{\displaystyle \frac{3}{s_{\mathrm{corr}}^2}}{\displaystyle _{L/2}^{L/2}}dz^{}\rho (z^{}){\displaystyle _0^{\mathrm{}}}d\varrho \varrho \left({\displaystyle \frac{1}{\sqrt{\varrho ^2+(zz^{})^2}}}{\displaystyle \frac{1}{s_{\mathrm{corr}}}}\right)\mathrm{\Theta }\left(s_{\mathrm{corr}}\sqrt{\varrho ^2+(zz^{})^2}\right)`$ (35)   
$`=`$ $`{\displaystyle \frac{3}{s_{\mathrm{corr}}^2}}{\displaystyle _{z_<}^{z_>}}dz^{}\rho (z^{}){\displaystyle _0^{\sqrt{s_{\mathrm{corr}}^2(zz^{})^2}}}d\varrho \varrho \left({\displaystyle \frac{1}{\sqrt{\varrho ^2+(zz^{})^2}}}{\displaystyle \frac{1}{s_{\mathrm{corr}}}}\right)`$ (36)   
$`=`$ $`{\displaystyle \frac{3}{2s_{\mathrm{corr}}^3}}{\displaystyle _{z_<}^{z_>}}dz^{}\rho (|z^{}|)\left(s_{\mathrm{corr}}|zz^{}|\right)^2,`$ (37)   
$``$

where $`z_<\mathrm{max}(L/2,zs_{\mathrm{corr}})`$, $`z_>\mathrm{min}(L/2,z+s_{\mathrm{corr}})`$ and $`s_{\mathrm{corr}}`$ is a function of $`z`$ through $`\rho (z)`$.

## IV Results and Conclusions

Once $`f_{\mathrm{corr}}`$, $`\mu _{\mathrm{ex}}(z)`$ and $`\rho _w(z)`$ are defined, the electric field and, consequently, the optimum density profile can be determined from the numerical iteration of Eq. (28) until convergence is obtained. The Helmholtz free energy, $`F`$, associated with the optimum counterion distribution (23), is determined by substituting it into the free-energy functional $``$,

$``$

$$\frac{\beta F}{𝒜}=\frac{2\sigma }{e}\left[\mathrm{ln}\left(\mathrm{\Lambda }^3\rho _0\right)1\right]\frac{\beta e}{2}_{L/2}^{L/2}dz\rho (z)\varphi (z)\beta \sigma \varphi \left(\frac{L}{2}\right)_{L/2}^{L/2}dz\rho (z)\left\{\beta \mu _{\mathrm{ex}}(z)\beta f_{\mathrm{corr}}\left[\rho _w(z)\right]\right\}.$$
(38)   
$``$

Using this expression, the pressure, for different distances between the plates, $`L`$, and various charge densities, $`\sigma `$, can be easily obtained through numerical differentiation,

$$P=\frac{1}{𝒜}\frac{F}{L},$$
(39)   
as shown in Fig. 2. When the charge density is below a threshold value, $`\overline{\sigma }<\overline{\sigma }_c`$, the dimensionless pressure, $`\lambda _B^3\beta P`$, is always positive and a monotonically decreasing function of $`\overline{L}`$. Above the critical surface charge density the pressure exhibits a distinct minimum. In particular, we find that for sufficiently high surface charge densities the force between the two like-charged surfaces becomes negative, i.e. the two plates attract!

In order to compare our results with other theories, we assumed that the dielectric medium between the plates is water at room temperature and, consequently, that the Bjerrum length is $`\lambda _B=7.14`$Å. The distance between the plates is fixed at 150 Å and the inverse of the surface charge density, $`\mathrm{\Sigma }=e/\sigma `$, is varied from 40 Å<sup>2</sup> to 1000 Å<sup>2</sup>. Our results, illustrated in Fig. 3, show that for small surface charge densities the pressure increases almost linearly with the inverse charge density, $`\mathrm{\Sigma }`$. In this case, since $`P_{\mathrm{corr}}P_{\mathrm{PB}}`$, the pressure is dominated by the PB behavior. However, when the charge density becomes large, the slope of $`P_{\mathrm{WDA}}`$ increases due to the strong repulsion between the counterions.

We also compare our calculations with the simulations of Guldbrand et al. In this case, the distance between the plates is fixed at 21 Å and the surface charge density is varied from 0.01 C/m<sup>2</sup> to 0.6 C/m<sup>2</sup>, as shown in Fig. 4.

When the density of counterions is small, $`P_{\mathrm{WDA}}`$ does not differ significantly from $`P_{\mathrm{PB}}`$. As the surface charge density is increased, the correlations among the counterions become relevant and $`P_{\mathrm{WDA}}`$ changes its slope and begins to decrease. Our results are in good agreement with the simulations, which also indicate that for a separation of 21 Å the pressure exhibits a region where it decreases with increase in the surface charge density.

###### Acknowledgements.

We acknowledge the fruitful discussions with Marcelo Louzada-Cassou, Rudi Podgornik, Roland Kjellander and Stjepan Marčelja. One of us, Marcia Barbosa, is particularly grateful for the useful discussion with Mark O. Robbins. This work was supported in part by CNPq — Conselho Nacional de Desenvolvimento Científico e Tecnológico and FINEP — Financiadora de Estudos e Projetos, Brazil. This research was also supported by the National Science Foundation under Grant No. PHY94-07194.