# Ion-ion correlations: an improved one-component plasma correction

## Abstract

Based on a Debye-Hückel approach to the one-component plasma we propose a new free energy for incorporating ionic correlations into Poisson-Boltzmann like theories. Its derivation employs the exclusion of the charged background in the vicinity of the central ion, thereby yielding a thermodynamically stable free energy density, applicable within a local density approximation. This is an improvement over the existing Debye-Hückel plus hole theory, which in this situation suffers from a “structuring catastrophe”. For the simple example of a strongly charged stiff rod surrounded by its counterions we demonstrate that the Poisson-Boltzmann free energy functional augmented by our new correction accounts for the correlations present in this system when compared to molecular dynamics simulations.

The classical one-component plasma (OCP) is an idealized model, in which a single species of ions moves in a homogeneous neutralizing background of opposite charge and interacts only via a repulsive Coulomb potential . Apart from its applications in plasma physics  it is also commonly used in soft matter physics as one of the simplest possible approaches for modeling correlations when studying polyelectrolytes, charged planes  or charged colloids . The general idea is the following: Compute the OCP free energy as a function of bulk density $`n_\mathrm{B}`$ and use this expression in the spirit of a local density approximation (LDA) as a correlation correction for the inhomogeneous system (i.e., $`n_\mathrm{B}n(𝒓)`$). The total excess free energy is the volume integral over the free energy density and thus becomes a functional of $`n(𝒓)`$. Many alternative and more sophisticated methods based on integral equations  have been developed for treating this correlation problem. Even though they offer results which are in good agreement with Monte-Carlo simulations, they do not provide any intuitive insight into the physics governing ionic solutions.

There is, however, a fundamental problem with the local density approaches: The OCP free energy is not a convex function of density . This implies that it cannot be used in a thermodynamically stable way within LDA, since the system can lower its total free energy by developing local inhomogeneities and increasing its density in one region at the expense of another (disregarding any surface effects) . Once started, this continues as a runaway process and the overall system collapses to a point. This feature is already seen on the level of the Debye-Hückel plus hole (DHH) approximation , which is an extension of the original Debye-Hückel (DH) theory  for the special case of the OCP, and the instability it gives rise to has been termed “structuring catastrophe” in this context . The proper way for avoiding this difficulty thus requires modifications of the one-component plasma model itself. The new theory, referred to as the Debye-Hückel-Hole-Cavity (DHHC) approach, remains simple and can be used within LDA to account for correlation effects present in more complex ionic solutions, as will be shown in an example at the end of the paper, where we compare its predictions to simulational results of a model system.

Since the necessary changes to DHH will turn out to be surprisingly tiny, it is worthwhile to briefly recall the way in which DHH theory arrives at a free energy for the OCP. For definiteness, we assume a system of $`N`$ identical point-particles of valence $`v`$ and (positive) unit charge $`q`$ inside a volume $`V`$ with a uniform neutralizing background of density $`vn_\mathrm{B}`$ and dielectric constant $`\epsilon `$. According to the DH approach, the potential $`\varphi `$ created by a central ion (i.e., fixed at the origin) and all its surrounding ions results from solving the spherically symmetric Poisson equation

$$^2\varphi (r)=\varphi ^{\prime \prime }(r)+\frac{2}{r}\varphi ^{}(r)=\frac{4\pi }{\epsilon }\rho (r)$$
(1)   
under the requirement that the charge density is $`\rho (𝒓)=qv\delta (𝒓)`$ at the central ion and that the rest of the mobile ions rearrange themselves in the uniform background in accordance with the Boltzmann distribution $`\rho (r)=vqn_\mathrm{B}\mathrm{exp}[\beta vq\varphi (r)]vqn_\mathrm{B}`$. Combining this with eq. (1) yields the nonlinear Poisson-Boltzmann (PB) equation, while linearization of the exponential function in the mobile ion density gives $`\rho (r)=\epsilon \kappa ^2\varphi (r)/4\pi `$ together with the famous Debye-Hückel solution for the potential, $`\varphi (r)\mathrm{e}^{\kappa r}/r`$, illustrating the rearrangement of the other ions around the central one in order to screen the Coulomb interaction. Here, $`\kappa \sqrt{4\pi \mathrm{}n_\mathrm{B}}`$ is the inverse screening length, $`\mathrm{}=\mathrm{}_\mathrm{B}v^2`$, with $`\mathrm{}_\mathrm{B}=\beta q^2/\epsilon `$ being the Bjerrum length, and $`\beta =1/k_\mathrm{B}T`$.

The problem with the DH theory is that the condition for linearization is obviously not satisfied for small $`r`$, where the potential is large — indeed, the particle density becomes negative and finally diverges at the origin. This defect was overcome by the DHH theory , which artificially postulates a correlation hole of radius $`h`$ around the central ion where no other ions are allowed. In this case the charge density is given by

$$\rho (r)=\{\begin{array}{cc}\hfill qv(\delta (𝒓)n_\mathrm{B}):& rh\hfill \\ \hfill \epsilon \kappa ^2\varphi (r)/4\pi :& r>h.\hfill \end{array}$$
(2)   
The solution of the linearized PB equation with the appropriate boundary conditions (continuity of electric field and potential) yields the potential for both regions in dependence of $`h`$, which has to be fixed on physical grounds: At low temperatures the electrostatic repulsion dominates and the minimum ion separation essentially becomes the mean separation, so $`h=(4\pi n_\mathrm{B}/3)^{1/3}`$. At high temperatures, the hole size can be estimated by balancing Coulombic and thermal energy, which gives $`h=\mathrm{}`$. A systematic way to interpolate between these two limits results from excluding particles from a region where their potential energy is larger than some threshold. A natural choice is the thermal energy $`k_\mathrm{B}T`$, which leads to

$$\kappa h=\omega 1\text{with}\omega =(1+3\mathrm{}\kappa )^{1/3}.$$
(3)   
Incidentally, this assumption also gives a continuous charge density across the hole boundary.

Once the potential at the position of the central ion is known, the electrostatic contribution to the Helmholtz free energy density can be obtained by the Debye charging process , as was done previously by Penfold et al. :

$$\frac{\beta f_{\mathrm{DHH}}}{n_\mathrm{B}}=\frac{1}{4}\left[1\omega ^2+\frac{2\pi }{3\sqrt{3}}+\mathrm{ln}\left(\frac{\omega ^2+\omega +1}{3}\right)\frac{2}{\sqrt{3}}\mathrm{arctan}\left(\frac{2\omega +1}{\sqrt{3}}\right)\right].$$
(4)   
The presented simple DHH analysis of the one-component plasma theory offers considerable insight into ionic systems and is in good agreement with Monte-Carlo simulations  when fluctuations on the charge density are not relevant . In principle one can attempt to include such fluctuations by applying the bulk density-functional theory in a local way, i.e., $`n_\mathrm{B}n(𝒓)`$. The basic idea is to obtain the density distribution via functional minimization of the Helmholtz free energy

$$\beta F_{\mathrm{OCP}}[n(𝒓)]=\mathrm{d}^3r\left\{n(𝒓)\mathrm{ln}\left(n(𝒓)V_\mathrm{p}\right)+\beta f_{\mathrm{DHH}}[n(𝒓)]\right\}$$
(5)   
under the constraint of global charge neutrality ($`V_\mathrm{p}`$ represents the volume of a particle.) Yet, this variational process does not lead to a well defined density profile, since $`f_{\mathrm{DHH}}(n)`$ asymptotically behaves like $`n^{4/3}`$ at high densities and therefore is not a convex function – with the implications already mentioned in the introduction. At small densities, however, the free energy density is convex and changes to a concave form only beyond a critical density $`n^{}7.8618/\mathrm{}^3`$ (see fig. 1). Hence, if during the process of actually computing $`n(𝒓)`$ such a density is never met, the theory does not “realize” its asymptotic instability and gives a finite (yet, meta-stable) answer. It has in fact been applied to account for correlations in the case of systems with low ionic strength ). Assuming the case of aqueous solutions ($`\mathrm{}_\mathrm{B}=7.14`$Å) and monovalent ions we find a critical density $`n^{}\mathrm{\hspace{0.33em}36}\text{mol/l}`$, which clearly is high enough to prevent a runaway process to set in. However, already for divalent and trivalent ions we find $`n^{}0.56\text{mol/l}`$ and $`0.049\text{mol/l}`$, respectively, which are sufficiently low to be realized and thus to trigger a collapse. Notice the strong dependence of $`n^{}`$ on valence, namely, on the sixth power.

To circumvent the instabilities occurring at high densities erroneously attributed to the local density approach itself, a number of nonlocal free energies have been proposed . In these weighted density approximations (WDA) the local density is replaced by a spatially averaged quantity. The main problem with these methods is that the choice of the weighting function is somewhat arbitrary. In most cases it is obtained by relating the second variation of the free energy with the direct correlation function. At this point the WDA requires prior information about this function, which is not yet available and thus has to be calculated using different approaches (like, e.g., integral equation theories). Whatever choice one takes, it is still (i) quite arbitrary and (ii) leads to a series of approximations which (iii) instead of clarifying the physics tend to obscure it.

The instabilities present in the local 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 LDA: 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. A natural solution for that problem is to decouple the particle density from the background density and to apply the LDA just to the former one. This, however, leads to nonlinearities in the solution of the differential equation which spoil the simplicity of the DH and DHH approximations. The most simple solution is to exclude the neutralizing background from a cavity of radius $`a`$ placed around the central ion only, which is already sufficient to control the unphysical divergence of the particle density. Even though it does not accounts for excluded-volume effects , the parameter $`a`$ can in principle be identified with the diameter of the particles. In addition, the exclusion hole for $`arh`$ is retained in order to account for the electrostatic repulsion between two ions. Consequently, the charge density, which for the usual DHH theory is given by eq. (2), has now three regions:

$$\rho (r)=\{\begin{array}{cc}\hfill qv\delta (𝒓):& r<a\hfill \\ \hfill qvn_\mathrm{B}:& arh\hfill \\ \hfill \epsilon \kappa ^2\varphi (r)/4\pi :& r>h.\hfill \end{array}$$
(6)   
The solution of the linearized PB equation with appropriated boundary conditions gives the potential in those regions:

$$\psi (r)=\frac{ve_0}{4\pi \epsilon r}\times \{\begin{array}{cc}1\frac{r}{2\mathrm{}}\left[(\kappa h)^2(\kappa a)^2\right]\kappa rC_h:\hfill & 0r<a\hfill \\ 1\frac{r}{2\mathrm{}}\left[(\kappa h)^2(\kappa r)^2\right]\frac{1}{3\mathrm{}\kappa }\left[(\kappa r)^3(\kappa a)^3\right]\kappa rC_h:\hfill & ar<h\hfill \\ C_h\mathrm{e}^{\kappa (rh)}:\hfill & hr<\mathrm{},\hfill \end{array}$$
(7)   
with the abbreviation

$$C_h=\frac{1}{1+\kappa h}\left(1\frac{(\kappa h)^3(\kappa a)^3}{3\mathrm{}\kappa }\right).$$
(8)   
In order to obtain the old theory in the limit $`a0`$ we choose the hole size $`h`$ to yield the same screening (i.e., the same amount of charge within $`h`$) as the DHH theory, which results in

$$\kappa h=\left[(\omega 1)^3+(\kappa a)^3\right]^{1/3}.$$
(9)   
This expression has four important physical limits: zero/infinite temperature and low/high density. At low temperature the exclusion hole has maximum size and, like in the DHH case, behaves as $`h=(3/4\pi n_\mathrm{B}+a^3)^{1/3}`$. As the temperature is increased, the hole size shrinks, but contrary to DHH theory it does not vanishes and $`ha`$ as $`T\mathrm{}`$. At small densities, entropic effects compete with the Coulombic repulsion and $`h=\mathrm{}+a`$; for high densities, the exclusion hole decreases but is again limited below and $`ha`$. Using this prescription for $`h`$, the Helmholtz free energy can be obtained by Debye-charging the fluid:

$$\frac{\beta f_{\mathrm{DHHC}}}{n_\mathrm{B}}=\frac{(\kappa a)^2}{4}_1^\omega \mathrm{d}\overline{\omega }\left\{\frac{\overline{\omega }^2}{2(\overline{\omega }^31)}\mathrm{\Omega }(\overline{\omega })^{2/3}+\frac{\overline{\omega }^3}{(1+\mathrm{\Omega }(\overline{\omega })^{1/3})(\overline{\omega }^2+\overline{\omega }+1)}\right\}$$
(10)   
with the abbreviation

$$\mathrm{\Omega }(\overline{\omega })=(\overline{\omega }1)^3+\frac{(\kappa a)^3}{3\mathrm{}\kappa }(\overline{\omega }^31)$$
(11)   
and where $`\omega `$ is the same as in eq. (3). The integral can be solved numerically for given values of $`\mathrm{}_\mathrm{B}`$, $`v`$ and $`a`$. As in the DHH approach fluctuations are taken into account by allowing the density to become local; thus, $`n(𝒓)`$ is obtained by minimizing the free energy from eq. (5) with $`f_{\mathrm{DHH}}`$ replaced by $`f_{\mathrm{DHHC}}`$ as given by eq. (10). But differently from the DHH theory, the Debye-Hückel-Hole-Cavity free energy is a convex function of density and thus applicable within a local density approximation. This situation is depicted in fig. 1, where we plotted the previous expression of the free energy of the DHH theory together with the improved expression of the DHHC approach. Recall that the DHH free energy has a point of inflection at a critical density $`n^{}7.8618/\mathrm{}^3`$, which makes it unstable at high densities – particularly for multivalent ionic correlations, as is demonstrated in the right part of fig. 1.

As an example, we apply this free energy as a correlation correction in the theoretical description of the screening of a charged rod, which is a simple model of biologically relevant stiff polyelectrolytes like DNA, actin filaments or microtubules. Much of the thermodynamic behavior of these molecules is determined by the distribution of the counterions around the polyion. As a model system we take a rod of radius $`r_0`$ and line charge density $`\lambda =0.959q/r_0`$ embedded in a cell of outer radius $`R=123.8r_0`$ and the complementary values $`\mathrm{}_\mathrm{B}/r_0=3`$, $`v=1`$ and $`\mathrm{}_\mathrm{B}/r_0=1`$, $`v=3`$ have been investigated, which on the plain PB level both give a fraction of condensed counterion (in the Manning sense) of roughly 65% . This system is thus strongly charged and one expects ionic correlations to become relevant. Indeed, the comparison between the distributions obtained by simulation  and the ones from PB theory shows that the mean-field approach fails in the limit of high ionic strength. In reality the ions do not just interact with the average electrostatic field but if an ion is present in a position $`𝒓`$, it tends to push away other ions from that point. This effect becomes important at high densities, low temperatures and for multivalent ions.

As discussed above, a simple way to improve PB theory is to extend the density functional to include a term of the form (10) which accounts for the correlations. The configurational free energy for the screened macroion solution can be partitioned into two terms:

$$F_\mathrm{P}[n(𝒓)]=F_{\mathrm{PB}}[n(𝒓)]+\mathrm{d}^3rf_{\mathrm{DHHC}}[n(𝒓)].$$
(12)   
The first part 

$$F_{\mathrm{PB}}[n(𝒓)]=\mathrm{d}^3r\left\{k_\mathrm{B}Tn(𝒓)\mathrm{ln}\left(n(𝒓)V_\mathrm{p}\right)+\frac{1}{2}qvn(𝒓)\varphi [n(𝒓)]\right\}$$
(13)   
contains the ideal gas contribution of the small ions, the interaction with the macroion potential and the mean-field interaction between the counterions. Minimization of this expression under the constraint of global charge neutrality gives – together with the Poisson equation – the Poisson-Boltzmann equation. The inter-particle correlations are now approximately accounted for by adding an excess free energy, which is the second term in eq. (12) — the DHHC free energy in local density approximation. The equilibrium ion distribution minimizing the functional (12) is most easily found by means of a Monte-Carlo solver, as has been proposed elsewhere . The fraction of ions within a distance $`r`$,

$$P(r)=\frac{1}{\lambda }_{r_0}^r\mathrm{d}\overline{r}\mathrm{\hspace{0.33em}2}\pi \overline{r}vqn(\overline{r}),$$
(14)   
obtained following this procedure is illustrated in fig. 2. Compared to the plain PB result the simulation shows a stronger condensation of ions in the vicinity of the rod, an effect which is more pronounced in the trivalent system. In both cases the increased condensation is reproduced by the correlation corrected PB functional from eq. (12). While in the case $`\mathrm{}_\mathrm{B}/r_0=3`$, $`v=1`$ the theoretical prediction practically overlaps the simulation, it somewhat overestimates correlations in the complementary case $`\mathrm{}_\mathrm{B}/r_0=1`$, $`v=3`$. It must, however, be noted that the ions in the simulation also interacted via a repulsive Lennard-Jones potential, giving them a diameter of roughly $`r_0`$. The expected reduction of particle density resulting from the additional hard core is not accounted for in the presented theory, but could easily be included along the lines of Refs. .

In conclusion, we have shown that the failure of the local density approximation for the one-component plasma is due to the asymptotically concave free energy employed by the DHH theory. To eliminate this problem, we introduced a DHHC approach in which the uniform background is absent in the immediate vicinity of the central ion, which leads to a convex, thermodynamically stable free energy. Moreover, the local density functional theory derived from this assumption is able to correctly account for the correlations between small ions in the presence of a strongly charged macroion. This was demonstrated for the case of a stiff rodlike polyelectrolyte by comparing the integrated charge density to simulation results of the same model. A more detailed investigation of the applicability of the LDA to rodlike polyelectrolytes and to charged colloids is postponed to future work.

###### Acknowledgements.

This work has been supported by the Brazilian agency CNPq (Conselho Nacional de Desenvolvimento Científico e Tecnológico). M. C. B. would like to thank K. Kremer for his hospitality during a stay in Mainz, where most of this work was completed.