Electrostatic correlations near charged planar surfaces

The fundamental description of charged systems is based on the well-known Poisson-Boltzmann (PB) theory, which was developed a century ago by Gouy¹ and Chapman.² 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_BT) 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⁷ 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²² within the Debye closure of the Bogoliubov-Born-Green-Kirkwood-Yvon hierarchical equations. Further, Netz and Orland²³ 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

$ü$ckel length κ⁻¹ = (8πq²l_Bρ_b)^−1/2 which represents the salt screening effects in the bulk; the Gouy-Chapman length l_G = 1/(2πl_Bqσ_s), describing the interaction between ion charge and surface charge; and the Bejerrum length l_B = e²/(4πε₀k_BT) representing the elementary charge interactions. Here, ρ_b, e, k_BT, and ε₀ 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²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,²⁴ 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.

In this section, we review the derivation of SC equations with the field theoretic variational approach.²³ 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 jth 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_BT), 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

$vc−1(r,r′)=−1βe2∇·(ε(r)∇δ(r−r′))$⁠, and the second part is self-energy of mobile ions, which should be subtracted from the total electrostatic energy with

$v_c^b(r) = l_B/r$

$vcb(r)=lB/r$ the bare Coulomb potential. Here, we should note that the self-energy term is crucial for electrostatic correlations. As a first step, the Lennard-Jones (LJ) interactions are ignored for point charge ions, we will consider the LJ interactions or other short range interactions in our finite size ions model in the future. The charge density operator is introduced as

$\rho ({\bf r}) = \sum _{i=1}^M \sum _{j=1}^{N_i} q_i \delta ({\bf r}-{\bf r}_{ij})+\rho _f({\bf r})$

$ρ(r)=∑i=1M∑j=1Niqiδ(r−rij)+ρf(r)$⁠, by assuming that all the ions are point-like charges. To account for the ion size effects into the system, a general charge distribution function, i.e., a Gaussian distribution function h(r, r^′), instead of δ(r, r^′), can be used for finite radius ions.⁸ The grand-canonical partition function

$\mathcal {Z}$

$Z$ is then transformed into a field-theoretic representation by the usual Hubbard-Stratonovich transformation, i.e.,

$\mathcal {Z}$

$Z$ is given by taking the form of a functional integral over the pure imaginary electrostatic auxiliary field ϕ(r); hence,

$\mathcal {Z} = \int \mathcal {D}\phi e^{-\mathcal {S}[\phi ]}$

$Z=∫Dϕe−S[ϕ]$⁠, with the action functional defined as

with

$λi=eiμ/λt3$ as the fugacity of ith 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

$\frac{\delta \mathcal {S}[\phi ]}{\delta \phi } = 0$

$δS[ϕ]δϕ=0$⁠, then it yields the MF theory equations. Furthermore, for a weakly correlated system, a perturbative expansion around the mean field gives more accurate results, however, it breaks down for strongly correlated systems. Another common way to evaluate the partition function or equivalent equilibrium measure

$dP(\phi )\break = \frac{1}{\mathcal {Z}} \exp (-\mathcal {S}(\phi ))$

$dP(ϕ)=1Zexp(−S(ϕ))$ is to create a stochastic process generated by the Langevin equation

$d\phi (\tau ) = -\frac{\partial \mathcal {S}}{\partial \phi (\tau )} d\tau + dW(\tau )$

$dϕ(τ)=−∂S∂ϕ(τ)dτ+dW(τ)$⁠, with 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

$\mathcal {S}_0$

$S0$⁠, the partition functional can be expressed as

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

$F=F0+⟨S−S0⟩0$⁠, where averages ⟨...⟩₀ are to be evaluated with respect to the most general Gaussian action function,

and with

$F0=−12trlnv0$⁠. The variational Gibbs free energy

$\mathcal {F}$

$F$ is then given by

By extremizing

$F$ with respect to v₀ and ϕ₀, i.e.,

$\delta \mathcal {F}/\delta v_0^{-1}({\bf r},{\bf r^{\prime }}) = 0$

$δF/δv0−1(r,r′)=0$ and

$\delta \mathcal {F}/\delta \phi _0({\bf r}) =0$

$δF/δϕ0(r)=0$⁠, we obtain the SC equations. For simplicity, only two ion species with symmetric valency are considered here:²³ 

where

$δw(r)=limr→r′[v0(r,r′)−vcb(r,r′)]$ and ε(r) = ε(r)/ε₀. We further rescale the SC equations with l_G, and by defining ψ = qϕ₀, ρ_f = ρ_fl_G/σ_s, v₀ = v₀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

$Λ=8πλlG3Ξexp(−ΞκlG/2)$⁠. Here, it is worth mentioning that we define the ionic self-energy dressed with electrostatic correlations in the form

$\delta v({\bf r}) = v_0({\bf r},{\bf r}) - v_c^b({\bf r},{\bf r})\break + l_B\kappa = v_0({\bf r},{\bf r}) - v_0^b({\bf r},{\bf r})$

$δv(r)=v0(r,r)−vcb(r,r)+lBκ=v0(r,r)−v0b(r,r)$⁠, with

$ v_0^b({\bf r}) = \frac{l_B}{r}\exp (-\kappa r)\break = \lim _{{\bf r} \rightarrow \infty } v_0({\bf r},{\bf r})$

$v0b(r)=lBrexp(−κr)=limr→∞v0(r,r)$⁠, the Debye-H

$\ddot{u}$

$ü$ckel potential in the homogeneous salt solution (i.e., Ξ = 0, ψ = 0).^3,4 That is, the reference system is chosen as the far field, such that we can obtain v₀(r, r) and

$v_0^b({\bf r})$

$v0b(r)$ simultaneously with same numerical accuracy, thus, the numerical errors for δv(r) are exactly canceled.

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

$ü$ckel equation that accounts for the local screening of the electrostatic propagator by mobile ions. Before we proceed to numerically solve the SC equations, we discuss their properties and subsequently we develop our numerical method based on these properties. First of all, we note that these two equations can only be decoupled in the limits Ξ → 0 and Ξ → ∞; and that an iterative method is needed for general Ξ. Second, these equations resemble the typical “Bratu problem” Δu + αe^u = 0,²⁶ 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 6D (r, r^′) modified Debye-H

$\ddot{u}$

$ü$ckle equation; however, only the diagonal information (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}$

$ü$ckel operator. Since the discrete Debye-H

$\ddot{u}$

$ü$ckel operator is a symmetric and positive denite matrix, we can use Cholesky or LDL factorization to fully invert the matrix, but the computational time cost (⁠

$\mathcal {O}(n^3)$

$O(n3)$⁠) and storage

$\mathcal {O}(n^2)$

$O(n2)$⁠) are extremely expensive for multi-dimensional problems. However, the problem is greatly simplified if only the diagonal elements are needed, which can be achieved by using the recently developed SelInv package to extract the diagonals of a matrix inverse,²⁴ with much smaller computational time cost (⁠

$\mathcal {O}(n^2)$

$O(n2)$⁠) and storage

$\mathcal {O}(n)$

$O(n)$⁠).

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

$ü$ckel operator for 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

$|\mathcal {F}^{(k)}-\mathcal {F}^{(k-1)}| < \delta$

$|F(k)−F(k−1)|<δ$⁠, with 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

$ü$ckel operator. The former is often approximated by the linearized PB equation, obtained by taking sinh (ψ(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,

$\Omega \equiv \cup _{j=1}^{M} \tau ^j$

$Ω≡∪j=1Mτj$ with τ^j the jth 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.0 ln (cosh (Δx)) + O((Δx)⁴) (Mickens discretization),²⁹ 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₀ the Bessel function. Then, the modified Debye-H

$ü$ckel equation can be decoupled as a series of one-dimensional equations for each k,

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

$|r‖−r‖′|=0$⁠, and the integral is obtained numerically via Gauss quadrature integration by a cutoff at a certain number, which is sufficiently large to account for the contribution from the smallest correlation length.

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

$ρ=∫0lG(ρ+(r)−ρ−(r))dr$⁠.

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.¹¹ 

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

$Π(h)=−∂F∂h$ as a function of separation distance. The results are shown in Figure 5. For a weakly charged system, the access free energy decreases by increasing the separation distance 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.

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)³³ 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,⁸ 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.