Electrostatic correlation effects near charged planar surfaces immersed in a symmetric electrolytes solution are systematically studied by numerically solving the nonlinear six-dimensional electrostatic self-consistent equations. We compare our numerical results with widely accepted mean-field (MF) theory results, and find that the MF theory remains quantitatively accurate only in weakly charged regimes, whereas in strongly charged regimes, the MF predictions deviate drastically due to the electrostatic correlation effects. We also observe a first-order like phase-transition corresponding to the counterion condensation phenomenon in strongly charged regimes, and compute the phase diagram numerically within a wide parameter range. Finally, we investigate the interactions between two likely-charged planar surfaces, which repulse each other as MF theory predicts in weakly charged regimes. However, our results show that they attract each other above a certain distance in strongly charged regimes due to significant electrostatic correlations.

## I. INTRODUCTION

The fundamental description of charged systems is based on the well-known Poisson-Boltzmann (PB) theory, which was developed a century ago by Gouy^{1} and Chapman.^{2} However, being a mean-field (MF) theory, the PB formalism ignores fluctuations and correlations, which are important for the cases of low temperature, highly charged surfaces, or multivalent counterions. In mesoscopic systems, the fluctuations are mainly due to the thermal noise, which is of the same order (*k*_{B}*T*) with other characteristic interactions. Also, the roughness of surfaces and solvent inhomogeneities etc. may also create some fluctuations. These fluctuation and correlation effects, which may drastically alter the mean-field picture of PB theory, have been the focus of recent theoretical efforts.^{3–10} For example, one surprising phenomenon beyond MF theory is the so-called fluctuation-driven counterion condensation in a simple system with single charged surface and its oppositely charged counterions.^{11–14} Another interesting phenomenon is the attraction between two highly likely-charged surface immersed in electrolytes or polar solution,^{14–18} as observed in experiments and in computer simulations. These phenomena cannot be captured within PB theories^{19–21} as PB theory only predicts pure repulsion between likely-charged surfaces. We will discuss these interesting phenomena further next.

A simple perturbative expansion about the MF solution near a charged surface by Netz and Orland^{7} demonstrated the breakdown of MF theory with relatively high charge densities. To go beyond the MF theory framework, nonlinear electrostatic self-consistent (SC) equations were first derived by Avdeev and Martynov^{22} within the Debye closure of the Bogoliubov-Born-Green-Kirkwood-Yvon hierarchical equations. Further, Netz and Orland^{23} re-derived the SC equations within the field theoretic formulation, using the Gaussian field variation ansatz. However, as stressed by the authors, these coupled nonlinear SC equations are computationally very complex. The main difficulties come from the high dimensionality (6D) and the nonlinearity of the coupled equations. Here, we propose a numerical method to solve these equations at an affordable computational cost based on the specific properties of the SC equations.

We consider relative simple systems with positive charged planar surfaces immersed in a symmetric electrolytes solution with fixed surface charge densities σ_{s} and salt valence *q*. The solvent is modeled implicitly as a continuous medium and the electric permittivity is a constant in the solution, hence, ignoring the solvent inhomogeneities in space, i.e., these inhomogeneity effects are considered as fluctuations. In general, there are three different length scales in such charged systems, i.e., the Debye-H

^{−1}= (8π

*q*

^{2}

*l*

_{B}ρ

_{b})

^{−1/2}which represents the salt screening effects in the bulk; the Gouy-Chapman length

*l*

_{G}= 1/(2π

*l*

_{B}

*q*σ

_{s}), describing the interaction between ion charge and surface charge; and the Bejerrum length

*l*

_{B}=

*e*

^{2}/(4πε

_{0}

*k*

_{B}

*T*) representing the elementary charge interactions. Here, ρ

_{b},

*e*,

*k*

_{B}

*T*, and ε

_{0}are the salt bulk concentration, elementary charge, thermal energy, and solvent dielectric permittivity, respectively. As suggested by previous studies,

^{14,23}we can define a key dimensionless parameter in terms of two lengths as Ξ =

*q*

^{2}

*l*

_{B}/

*l*

_{G}. On the one hand, the limit Ξ → 0 is called weak coupling limit, where the thermal fluctuations overcome electrostatic interactions, thus the physical system can be precisely described within MF or PB theory. On the other hand, Ξ → ∞ is the strong coupling limit, where the electrostatic correlations dominate the thermal fluctuations.

^{3,23}A perturbative approach around the MF theory may not be valid, instead, non-perturbative approaches (e.g., the variational approach) remain accurate in this regime.

The central focus of the current work is to provide accurate numerical results from solving the nonlinear electrostatic SC equations. The simulation data cover a wide parameter range, and hence we can examine a complete counterion condensation phase diagram. We further explore other important features of the counterion condensation transition, such as the order of the phase transition, as well as the scaling properties of the order parameter and free energy near the transition point. We also investigate the interactions between two likely-charged surfaces as a function of separation distance in salt solution with different coupling parameters.

The paper is organized as follows: we start by re-deriving of the SC equations via the field theory approach. We then solve the equations numerically (see Sec. III) with the help of selective inverse techniques,^{24} saving both memory storage and computational time. In Sec. IV, we describe the main numerical results obtained from the SC equations. In Sec. IV A, we describe the counterion condensation phase diagram and analyze the associated scaling behavior near the transition point. In Sec. IV B, we show the full interaction pictures between two-likely charged surfaces with electrostatic correlation effects. We conclude with a short discussion about the limitations and future studies are included in Sec. V.

## II. SELF-CONSISTENT (SC) EQUATIONS

In this section, we review the derivation of SC equations with the field theoretic variational approach.^{23} The grand canonical partition function of *M* ion species in a liquid of spatially varying dielectric permittivity ε(**r**) is

where **r**_{ij} is the position of *j*th ion of type *i*, λ_{t} is the thermal wavelength of point-like ions, and *N*_{i} and μ_{i} are the total number and chemical potential of ions of type *i*, respectively. The indicator function χ(**r**) takes into account the presence of hard walls, χ(**r**) = 1 in the solution, and χ(**r**) = 0 inside the hard walls. In general, the hard surface can be represented by a hard potential as χ(**r**) = exp ( − *V*(**r**)/*k*_{B}*T*), *V*(**r**) = ∞, and we can adjust *V*(**r**) to represent different surfaces. *H* is the Hamiltonian of the charged system, which consists of two parts as follows:

where the first part is the electrostatic energy with Coulomb potential defined as

*h*(

**r**,

**r**

^{′}), instead of δ(

**r**,

**r**

^{′}), can be used for finite radius ions.

^{8}The grand-canonical partition function

**r**); hence,

with

*i*th ion specie, here, we assume λ

_{+}= λ

_{−}= λ for symmetric salt solution. However, it is still impossible to calculate the above functional integral to obtain the partition function directly. Various approximations have been introduced in statistical physics to solve this problem. The simplest yet most popular approximation is to replace the functional integral over the electrostatic fields by the value of the integrand at the extremum or optimal value, i.e., the saddle point approximation, which determines the optimal values of the field given by

*W*(τ) the Wiener process, and with the ergodic stochastic process having the same measure

*dP*(ϕ) as

*t*→ ∞. However, this suffers from the long-time integration problem in solving numerically this stochastic partial differential equations. A more powerful and effective non-perturbative way is the standard Gibbs variational procedure. Based on the chosen proper reference action function

The above inequality holds due to Jensen's Inequality.^{25} Then, we approximate the partition function by the upper bound or optimizing the variational grand canonical free energy

_{0}are to be evaluated with respect to the most general Gaussian action function,

and with

By extremizing

*v*

_{0}and ϕ

_{0}, i.e.,

^{23}

where

*r*) = ε(

*r*)/ε

_{0}. We further rescale the SC equations with

*l*

_{G}, and by defining ψ =

*q*ϕ

_{0}, ρ

_{f}= ρ

_{f}

*l*

_{G}/σ

_{s},

*v*

_{0}=

*v*

_{0}

*l*

_{G}/

*l*

_{B}, we obtain the dimensionless equations as

Clearly, there are only two dimensionless parameters in this physical system from Eq. (8), namely, the coupling parameter Ξ defined in Sec. I, and the rescaled fugacity

^{3,4}That is, the reference system is chosen as the far field, such that we can obtain

*v*

_{0}(

**r**,

**r**) and

*v*(

**r**) are exactly canceled.

## III. NUMERICAL METHODS

The first of the SC equations is a modified PB equation for the fluctuating external potential induced by the fixed surface charge, while the second is a modified Debye-H

*u*+ α

*e*

^{u}= 0,

^{26}the highly nonlinear term

*e*

^{u}becomes a big trouble near bifurcation points for numerical methods. We note that most of the iterative methods, i.e., fixed point iteration method or Newton iteration method only work well before the bifurcation point. However, the bifurcation point is isolated, thus we can approach the point from one side of parameter space with the continuation method. Physically, the bifurcation point corresponds to the famous counterion condensation phase transition. Another important aspect and also the greatest difficulty in solving the SC equations arise from the 6

*D*(

**r**,

**r**

^{′}) modified Debye-H

**r**=

**r**

^{′}) is needed for calculating δ

*v*(

**r**). Thus,

*the problem reduces to finding the diagonal elements of the inverse of the discrete Debye-H$\ddot{u}$$u\u0308$ckel operator*. Since the discrete Debye-H

^{24}with much smaller computational time cost (

We employ a fixed-point iteration scheme for the solution of coupled SC equations following two alternating steps: first, for given δ*v*(**r**), we solve the modified PB equation for ψ(**r**) given the boundary conditions ψ(**r** → ∞) = 0; next for given δ*v*(**r**) and ψ(**r**), we invert the discrete Debye-H

*v*(

**r**,

**r**′), and then obtain a new δ

*v*(

**r**). These two steps are iteratively performed until the convergence criterion of the solution is reached; here, we use the relative free energy as the convergence criterion

*k*= 1, 2, … the iteration step and δ a small pre-set value.

We consider numerical solutions of the nonlinear modified PB equation and the discrete form of the Debye-H

**r**)) ≈ ψ(

**r**) when ψ(

**r**) ≪1. The nonlinear modified PB equation presents several numerical difficulties due to the rapid exponential nonlinearities, discontinuous coefficients, delta functions (point charge distribution), etc. The finite volume method has been one of the standard approaches for discretizing 2D and 3D interface problems which enforce conservation of certain physical quantities even for the discrete systems. First, we partition the domain Ω into 3D non-uniform Cartesian mesh,

^{j}the

*j*th finite volume. It is important to notice that all the discontinuities in the coefficients are taken to lie along the boundaries of τ

^{j}. Then, by integrating the modified PB equation over an arbitrary finite volume, and employing the divergence theorem for the first term in the equation, we obtain the algebraic equations in the form

**A**ψ +

*N*(ψ) −

*f*= 0. We then solve the algebraic equation with inexact-Newton method for the modified PB equation.

^{27,28}A small trick is used during discretization for the “Bratu” type problem, i.e., we approximate the second-order derivative using a nonlinear denominator function, (Δ

*x*)

^{2}= 2.0 ln (cosh (Δ

*x*)) +

*O*((Δ

*x*)

^{4}) (Mickens discretization),

^{29}which works well near bifurcation points.

For simple systems with charged planar interfaces located at the (*y*, *z*) plan, the translational symmetry will greatly simplify the problem. In such cases, the electrostatic potential becomes simply a function of the coordinate *x*, and the electrostatic Green's function can be expanded in 2D Fourier basis, with the Fourier transform of circularly symmetrical function (Hankel transform) as

with *J*_{0} the Bessel function. Then, the modified Debye-H

*k*,

Here δ*v*(*x*) is calculated using Eq. (10) with

Once we solve the above SC equations, the local ion densities are given by

To further investigate the physical phenomena near the charged surface, we define a proper order parameter ρ for ion adsorption by the charged surface as

## IV. RESULTS AND DISCUSSION

### A. Single planar surface in symmetric electrolytes solution

First, we consider one planar charged surface case with fixed charge density ρ_{f}(*x*) = δ(*x*). The results from numerically solving the SC equations with different coupling parameters, such as electrostatic potential, correlation potential, and ion density distribution near one charged surface, are shown in Figure 1. Clearly, when Ξ is small, the results from the SC equations approach the same results as the one from the PB equation. However, as Ξ increases, the results differ from the PB results, especially for larger value of Ξ. One of the interesting phenomena is that there is a peak in the positive ion density profile near the positive charged surface. This is due to the strong correlations, which cause overcompensation of surface charge by negative ions close to or at the surface. Then, the positive ions which are away from the surface will “feel” a negative apparent charged surface. Therefore, they will be attracted toward the wall. This is the so called “charge inversion” which was observed by previous studies.^{30–32}

As we already discussed at the end of Sec. II, the state of the simple charged systems is determined by two parameters: *the coupling parameter Ξ*, and *the rescaled fugacity Λ*. For various rescaled fugacity values, we change Ξ from small to large values with small increments and we measure the free-energy and the adsorbed amount of ions within a Gouy-Chapman length near the wall ρ. There is obviously a phase transition when Ξ become bigger corresponding to the bifurcation point of the SC equation. Above this point, simple fixed point iteration method will blow up. All of the data near the transition point are then carefully analyzed.

The numerical results for the phase transition between the weakly ions adsorption and counterion condensation states based on the analysis of the free energy are displayed in Figure 2. The phase boundary curve is constructed by simulation data within three orders of magnitude in parameter space. The fact that the counterion condensation transition is a first-order phase transition is established. As a result, going from the weakly adsorption state, the free energy and ion adsorbed amount follow some scaling behavior near the transition point. To demonstrate the first-order phase transition, we plot (ρ(Ξ, Λ) − ρ_{c}) as a function of (Ξ/Ξ_{c} − 1) in a double logarithmic plot in Figure 3(a). A dashed line representing slope of 1 has also been included. All data collapse and approach a single scaling power law as shown in Eq. (12):

To confirm the first-order nature of the transition, we further examine the dependence of the free energy per unit volume as a function of (Ξ/Ξ_{c} − 1) near the transition point. For a first-order transition, we expect that

The power law exponent 1 is characteristic of the first-order transition. Figure 3(b) demonstrates that all Λ cases display the above power law near the transition point; for comparison, a dashed line has been drawn with the anticipated power-law exponent 1. Our numerical results are consistent with previous studies.^{11}

### B. Two planar surfaces in symmetric electrolyte solution

We now consider the case of two charged planar surfaces immersed in salt solution with separation distance *h*, and fixed charged density ρ_{f}(*x*) = δ(*x* − 0.5*h*) + δ(*x* + 0.5*h*). First, we show electrostatic potential, correlation potential, and ion density profiles in Figure 4. For small separation distance *h*, or strong confinement, the osmotic pressure is extremely big, pushing the ions into the middle of the confinement layer, see black lines of Figures 4(g) and 4(h). The osmotic pressure is a key parameter in understanding the thermodynamics and the interactions in colloidal systems. Here, we calculate the osmotic pressure from

*h*, while the osmotic pressure is pure repulsion. However, in strongly charged systems, the access free energy decreases when we increase the separation distance only within a small regime, i.e.,

*h*∼

*l*

_{G}, and the access free energy then increases by increasing the separation distance for large

*h*, which corresponds to an attractive osmotic pressure as shown in Figure 5(d). Finally, by calculating all the osmotic pressure data in the parameter space, we obtain the global phase diagram for the interactions between two likely-charged surfaces as shown in Figure 6. The phase boundary shape here is similar to the one of the counterion condensation phase diagram (Figure 2) with a shift on the log-log plot due to the confinement effects.

## V. CONCLUSIONS

In this work, we have developed a numerical method for solving the electrostatic SC equations in order to explore the parameter space for simple charged planar surfaces immersed in salt solution. We quantified the physical picture beyond MF theory, i.e., we established the first-order phase transition properties of counterion condensation and the phase diagram, and we calculated the osmotic pressure between two likely-charged planar surface. For the latter case, we found that the surfaces may attract each other for a certain separation distance in strongly charged systems.

For the moment, we only consider very simple geometry systems, however, the extension to complex geometry systems is straightforward. From our current studies, we observe that the electrostatic correlation energies are mainly located near the charged surface, thus, an adaptive mesh refinement scheme should be constructed with most of the grid points distributed near the charged surfaces. Alternatively, a high-order Discontinuous Galerkin (DG) method or extensions of the finite volume method employed here may result in higher accuracy and efficiency.

Due to the bifurcation nature of the SC equations, our current fixed point iteration method will blow up above the bifurcation point, thus, there remains one problem, namely, what should we do after the bifurcation point? To this end, we propose the recursive projection method (RPM)^{33} to deal with these situations. On the other hand, physically, we can always renormalize the charge densities and the condensation ions into a new charge density, such that the renormalized system can be solved with our current method again.

Another important aspect is to account for the ion size effects^{34–36} in SC equations,^{8} i.e., the charge on an ion is assumed to have a finite spread around its center, instead of point charge δ(**r** − **r**^{′}). A convenient choice is the Gaussian distribution with Born radius *a*,

The LJ interactions are ignored for point charges here as a first step, we plan to include the LJ interactions or other short range interactions in our finite size ions model. We will discuss the new model and corresponding numerical results in our future work.

## ACKNOWLEDGMENTS

This work was sponsored by the Collaboratory on Mathematics for Mesoscopic Modeling of Materials (CM4) supported by DOE, and also by NIH Grant No. 1U01HL114476-01A1. Computations were performed using a DOE INCITE award. We would like to acknowledge Dr. Zhongqiang Zhang of Brown University for helpful discussions.