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.
The fundamental description of charged systems is based on the well-known Poisson-Boltzmann (PB) theory, which was developed a century ago by Gouy1 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 (kBT) 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 theories19–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 Orland7 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 Martynov22 within the Debye closure of the Bogoliubov-Born-Green-Kirkwood-Yvon hierarchical equations. Further, Netz and Orland23 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
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 rij is the position of jth ion of type i, λt is the thermal wavelength of point-like ions, and Ni 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)/kBT), 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
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
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
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
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
We consider numerical solutions of the nonlinear modified PB equation and the discrete form of the Debye-H
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 J0 the Bessel function. Then, the modified Debye-H
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.5h) + δ(x + 0.5h). 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
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 effects34–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.
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.