Mao Su(苏茂)1,2, Yanting Wang(王延颋),1,21CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, 55 East Zhongguancun Road, PO Box 2735, Beijing 100190, China 2School of Physical Sciences, University of Chinese Academy of Sciences, 19A Yuquan Road, Beijing 100049, China
Abstract The Poisson–Boltzmann (PB) theory is one of the most important theoretical models describing charged systems continuously. However, it suffers from neglecting ion correlations, which hinders its applicability to more general charged systems other than extremely dilute ones. Therefore, some modified versions of the PB theory are developed to effectively include ion correlations. Focused on their applications to ionic solutions, the original PB theory and its variances, including the field-theoretic approach, the correlation-enhanced PB model, the Outhwaite–Bhuiyan modified PB theory and the mean field theories, are briefly reviewed in this paper with the diagnosis of their advantages and limitations. Keywords：Poisson–Boltzmann theory;ion correlation;ionic solution;field-theoretic approach;mean field theory
PDF (484KB)MetadataMetricsRelated articlesExportEndNote|Ris|BibtexFavorite Cite this article Mao Su(苏茂), Yanting Wang(王延颋). A brief review of continuous models for ionic solutions: the Poisson–Boltzmann and related theories. Communications in Theoretical Physics, 2020, 72(6): 067601- doi:10.1088/1572-9494/ab8a23
1. Introduction
Theoretical descriptions of charged systems are very important for a wide range of applications, such as the interfacial phenomena of aqueous solutions [1–3] and the structure and dynamics of biomolecules [4–6]. Among various kinds of theories, the Poisson–Boltzmann (PB) theory initially developed by Debye and Hückel [7] in 1923 is the most famous one. The PB theory is particularly useful in depicting ionic solutions, for which it serves as a mesoscopic continuous model describing the average potential distribution of ionic solutions. Once the potential distribution is determined, physical properties such as activity and osmotic coefficient can be calculated with the knowledge of statistical physics. Because of their simplicity, the PB theory and its linearized version, the Debye–Hückel (DH) theory, are widely used in many scientific occasions, including investigating the electrostatic and thermodynamic properties, such as activity coefficient, solvation energy, and so on [8, 9], of ionic solutions and charged molecules.
It is well known that the PB theory is on the basis of key approximations we will introduce below that limit its applicability to extremely dilute solutions whose ion correlation effect is negligible. Therefore, the PB model can completely fail in some practical scenarios [10]. For example, DNA and many protein molecules are so highly charged that the electrostatic correlation effect cannot be neglected [6, 11]. It is worth noting that the PB theory may in some occasions, e.g. calculating thermodynamic properties for densely charged systems, provide reasonable results due to error cancellation. However, it is by chance and cannot serve as a routine treatment to dense systems.
As illustrated in figure 1, in order to yield more accurate solutions for a wider range of charged systems, a variety of theoretical models, either based on the PB theory or not, have been developed [12–14]. Interestingly, most of them are connected to the PB theory, either through the form of equation or through the form of solution, which in turn manifests the importance of the PB theory.
Figure 1.
New window|Download| PPT slide Figure 1.Schematic of a charged system which can be described by the PB theory. When the correlation effect illustrated by the halos around point charges is significant, the original PB theory may fail and enhanced PB theories are on demand.
This review is organized as follows. The derivation of the original PB theory is shown in section 2, along with detailed discussion of the approximations used in the derivation. In section 3, the field-theoretic (FT) approach for calculating the partition function is reviewed, and as a consequence a set of self-consistent equations are obtained by the variational method. In section 4, we introduce our correlation-enhanced PB model along with the Outhwaite–Bhuiyan modified PB (MPB) theories. Finally, a set of mean field theories including the dressed-ion theory (DIT) and the molecular Debye–Hückel (MDH) theory are reviewed in section 5. In the appendix, the integral equation theory is briefly introduced as the concepts are used in the MPB theories as well as the mean field theories. The SI unit is used by default unless the use of the Gaussian unit is explicitly declared. It should be noted that this brief review is more titled to serve as a manual for beginners, so some most recent advanced theoretical and numerical progresses, such as the concentration-dependent dielectric model [15, 16], are not reviewed due to the limited space.
2. The PB theory
The PB theory is originally derived by Debye and Hückel in 1923 [7]. It is formally a Poisson’s equation with the charge density distribution obeying the Boltzmann distribution. The original PB equation is a nonlinear equation which is usually hard to solve analytically. Debye and Hückel then implement some approximations to linearize the original PB equation, leading to the linearized PB equation, or the DH equation [17], which can be solved analytically and end up with the Yukawa potential. The activity coefficients calculated by the DH equation at the low-concentration limit fit the experimental values very well, but deviate more severely at a higher concentration due to the approximations incorporated in the DH equation [17].
2.1. Derivation of the original PB theory
The derivation of the PB theory can be found in [17] and we summarize it as follows. The first approximation made in the PB theory is neglecting the inner structure of particles. For an aqueous ionic solution, water is treated as a continuous background with the relative dielectric constant $\varepsilon ,$ and ions are assumed to be point charges. For an infinitely large domain consisted of N ions located at ${r}_{1},{r}_{2},\cdots ,{r}_{N},$ the electric potential of the system at $r$ is$\begin{eqnarray}\psi (r)=\displaystyle \sum _{i=1}^{N}\displaystyle \frac{{q}_{i}}{4\pi {\varepsilon }_{0}\varepsilon \left|r-{r}_{i}\right|},\end{eqnarray}$which satisfies Poisson’s equation$\begin{eqnarray}-{{\rm{\nabla }}}^{2}\psi (r)=\displaystyle \frac{1}{{\varepsilon }_{0}\varepsilon }\rho (r),\end{eqnarray}$where ${\varepsilon }_{{\rm{0}}}$ is the dielectric constant of vacuum, ${q}_{i}$ is the charge of the ith ion, and $\rho \left(r\right)$ is the charge density at r. In a homogeneous charge-neutral system, the solution of $\psi \left(r\right)$ is zero everywhere. To obtain useful information, we fix an ion at ${r}_{1}$ and study the electrostatic potential at $r:$$\begin{eqnarray}\psi (r,{r}_{1})=\displaystyle \frac{{q}_{1}}{4\pi {\varepsilon }_{0}\varepsilon \left|r-{r}_{1}\right|}+\displaystyle \sum _{i=3}^{N}\displaystyle \frac{{q}_{i}}{4\pi {\varepsilon }_{0}\varepsilon \left|r-{r}_{i}\right|}.\end{eqnarray}$
Since the positions of other ions fluctuate, we find the canonical ensemble average of $\psi (r,{r}_{1})$ to be:$\begin{eqnarray}\left\langle \psi (r,{r}_{1})\right\rangle =\displaystyle \frac{\displaystyle \int \cdots \displaystyle \int \psi (r\,,{r}_{1}){{\rm{e}}}^{-\beta U}{\rm{d}}{r}_{2}\cdots {\rm{d}}{r}_{N}}{\displaystyle \int \cdots \displaystyle \int {{\rm{e}}}^{-\beta U}{\rm{d}}{r}_{2}\cdots {\rm{d}}{r}_{N}},\end{eqnarray}$where $U=U\left({r}_{1},{r}_{2},\mathrm{...},{r}_{N}\right)$ is the potential energy and $\beta \equiv \tfrac{1}{{k}_{{\rm{B}}}T}$ with ${k}_{{\rm{B}}}$ the Boltzmann constant and $T$ the temperature. Taking the Laplacian with respect to $r$ on both sides:$\begin{eqnarray}\begin{array}{lll}{{\rm{\nabla }}}^{2}\left\langle \psi (r,{r}_{1})\right\rangle & = & \displaystyle \frac{\displaystyle \int \cdots \displaystyle \int {{\rm{\nabla }}}^{2}\psi (r,{r}_{1}){{\rm{e}}}^{-\beta U}{\rm{d}}{r}_{2}\cdots {\rm{d}}{r}_{N}}{\displaystyle \int \cdots \displaystyle \int {{\rm{e}}}^{-\beta U}{\rm{d}}{r}_{2}\cdots {\rm{d}}{r}_{N}}\\ & = & \displaystyle \frac{\displaystyle \int \cdots \displaystyle \int -\displaystyle \frac{{\rm{1}}}{{\varepsilon }_{0}\varepsilon }\rho (r,{r}_{1}){{\rm{e}}}^{-\beta U}{\rm{d}}{r}_{2}\cdots {\rm{d}}{r}_{N}}{\displaystyle \int \cdots \displaystyle \int {{\rm{e}}}^{-\beta U}{\rm{d}}{r}_{2}\cdots {\rm{d}}{r}_{N}}\\ & = & -\,\displaystyle \frac{1}{{\varepsilon }_{0}\varepsilon }\left\langle \rho (r,{r}_{1})\right\rangle .\end{array}\end{eqnarray}$
The average density distribution $\left\langle \rho (r,{r}_{1})\right\rangle $ is related to the radial distribution function (RDF) $g\left(r,{r}_{1}\right)$ by$\begin{eqnarray}\left\langle \rho (r,{r}_{1})\right\rangle =\displaystyle \displaystyle \sum _{s}{c}_{s}{q}_{s}{g}_{1s}(r,{r}_{1}),\end{eqnarray}$where ${c}_{s}$ is the bulk number concentration of the s-type ions and ${g}_{1s}\left(r,{r}_{1}\right)$ is the RDF of s-type ions about the central ion located at ${r}_{1}.$ Note that the mean potential $\left\langle \psi (r,{r}_{1})\right\rangle $ and the RDF ${g}_{1s}\left(r,{r}_{1}\right)$ depend only upon $\left|r-{r}_{1}\right|,$ we can use $r$ instead of $(r,{r}_{1})$ to simplify the notations.
The potential of mean force (PMF) was further utilized for finding the solution of the above equation [18]. We Consider the mean force exerted on the ion at ${r}_{2}$ with the fixed ion at ${r}_{1}:$$\begin{eqnarray}\begin{array}{lll}-\left\langle \displaystyle \frac{{\rm{d}}}{{\rm{d}}{r}_{2}}U\right\rangle & = & \displaystyle \frac{-\displaystyle \int \cdots \displaystyle \int \displaystyle \frac{{\rm{d}}U}{{\rm{d}}{r}_{2}}{{\rm{e}}}^{-\beta U}{\rm{d}}{r}_{3}\cdots {\rm{d}}{r}_{N}}{\displaystyle \int \cdots \displaystyle \int {{\rm{e}}}^{-\beta U}{\rm{d}}{r}_{3}\cdots {\rm{d}}{r}_{N}}\\ & = & \displaystyle \frac{1}{\beta }\displaystyle \frac{\displaystyle \int \cdots \displaystyle \int \displaystyle \frac{{\rm{d}}}{{\rm{d}}{r}_{2}}{{\rm{e}}}^{-\beta U}{\rm{d}}{r}_{3}\cdots {\rm{d}}{r}_{N}}{\displaystyle \int \cdots \displaystyle \int {{\rm{e}}}^{-\beta U}{\rm{d}}{r}_{3}\cdots {\rm{d}}{r}_{N}}\\ & = & \displaystyle \frac{1}{\beta }\displaystyle \frac{{\rm{d}}}{{\rm{d}}{r}_{2}}\,\mathrm{ln}\left(\displaystyle \int \cdots \displaystyle \int {{\rm{e}}}^{-\beta U}{\rm{d}}{r}_{3}\cdots {\rm{d}}{r}_{N}\right)\\ & = & \displaystyle \frac{1}{\beta }\displaystyle \frac{{\rm{d}}}{{\rm{d}}{r}_{2}}\,\mathrm{ln}\left(\displaystyle \frac{N(N-1)\displaystyle \int \cdots \displaystyle \int {{\rm{e}}}^{-\beta U}{\rm{d}}{r}_{3}\cdots {\rm{d}}{r}_{N}}{Z}\right)\\ & \equiv & \displaystyle \frac{1}{\beta }\displaystyle \frac{{\rm{d}}}{{\rm{d}}{r}_{2}}\,\mathrm{ln}\,g({r}_{2},{r}_{1}).\end{array}\end{eqnarray}$
In the above derivation, the ion number N and the partition function $Z=\displaystyle \int \cdots \displaystyle \int {{\rm{e}}}^{-\beta U}{\rm{d}}{r}_{1}\cdots {\rm{d}}{r}_{N}$ are independent of ${r}_{2},$ so $N\left(N-1\right)/Z$ can be inserted to obtain the definition of the RDF. We then define$\begin{eqnarray}w({r}_{2},{r}_{1})\equiv -\displaystyle \frac{1}{\beta }\,\mathrm{ln}\,g({r}_{2},{r}_{1}).\end{eqnarray}$
Obviously the negative gradient of $w({r}_{2},{r}_{1})$ is the mean force, so $w({r}_{2},{r}_{1})$ is called the PMF. Similarly, the PMF depends only upon the distance $\left|{r}_{2}-{r}_{1}\right|$, then we have $\mathrm{In}\,g(r)=-\beta w(r)$.
The second approximation is that the PMF ${w}_{1s}\left(r\right)$ is replaced by ${q}_{s}\psi \left(r\right),$ where the mean potential energy $\psi \left(r\right)\equiv \left\langle \psi \left(r\right)\right\rangle ,$ which leads to the PB equation:$\begin{eqnarray}{{\rm{\nabla }}}^{2}\psi (r)=-\displaystyle \frac{1}{{\varepsilon }_{0}\varepsilon }\displaystyle \displaystyle \sum _{s}{c}_{s}{q}_{s}{{\rm{e}}}^{-\beta {q}_{s}\psi (r)}.\end{eqnarray}$
When the third approximation of linearization (expanding the exponent to the first order) is taken, the PB theory is further simplified to be the analytically solvable DH equation:$\begin{eqnarray}{{\rm{\nabla }}}^{2}\psi (r)=\displaystyle \frac{\beta }{{\varepsilon }_{0}\varepsilon }\displaystyle \displaystyle \sum _{s}{c}_{s}{q}_{s}^{2}\psi (r).\end{eqnarray}$
The solution of the DH equation for a point charge in an infinitely large domain is a screened Coulomb potential, named the Yukawa potential:$\begin{eqnarray}\psi (r)=\displaystyle \frac{{q}_{1}}{4\pi {\varepsilon }_{0}\varepsilon r}{{\rm{e}}}^{-\kappa r},\end{eqnarray}$where the Debye screening length $1/\kappa $ is defined by$\begin{eqnarray}{\kappa }^{2}\equiv \displaystyle \frac{\beta }{{\varepsilon }_{0}\varepsilon }\displaystyle \displaystyle \sum _{s}{c}_{s}{q}_{s}^{2}.\end{eqnarray}$
The physical properties such as RDF, activity coefficient, and osmotic coefficient can be calculated once the potential is known [17].
The PB equation is difficult to be solved due to the nonlinearity. However, in some special cases it is possible to find the analytical solutions. For example, the PB theory can be applied to electrolytes with a charged plate, and the result is known as the Gouy–Chapman theory [19]. More generally, the PB equation is solved by numerical methods such as the finite difference method, finite element method, and so on. Accelerate techniques such as conjugate gradient method, fast multi-pole method, and fast Fourier transform-based approaches are also used to solve the PB equation numerically [20].
2.2. Limitations
The PB theory neglects molecular structures (exclusive volume effects) of ions and solvents. As a result, the structural properties calculated by the PB theory may deviate significantly from the real solutions. For example, the RDF for an ionic aqueous solution described by the PB theory is essentially the density distribution with respect to a hard sphere and the central ion is treated as a boundary condition, so the RDF has the feature of a gas, i.e. with only one peak and decaying monotonically [21]. In real solutions, however, a typical RDF always has a valley at about 0.5 nm, which is about the size of a water molecule, because that space is occupied explicitly by water molecules [22].
In most application scenarios, however, the exclusive volume effect only has a minor influence because thermodynamic properties such as solvation energy are the major focus of interest. Implicit solvent models are studied by researchers to evaluate the solute PMF, which determines the statistical weight of solute conformations and therefore allows molecular simulations without explicit solvent molecules [23]. When the solution goes to the high concentration limit, i.e. near close packing, the PB theory completely fails due to the significant exclusive volume effect. In this case, the integral equation theory (see appendix) with the Percus–Yevick (PY) or Hypernetted Chain (HNC) approximation can be employed to calculate the distribution function [24]. In particular, the HNC approximation is known to be very accurate for the charged hard-sphere model. However, the integral equation method cannot be applied to problems involving boundaries.
The second approximation of replacing the PMF by the potential energy $q\psi \left(r\right)$ in the PB theory neglects ion correlations, which is however most of the time important due to the long-range nature of electrostatic interaction. Many phenomena such as overcharge [25–29], like-charge attraction [30–34], and opposite-charge repulsion [35] have correlations and thus cannot be satisfactorily described by the PB theory. Particularly, the charge effect plays an essential role in most biomolecular systems, but many of them, such as protein and DNA [4–6, 25], cannot be effectively described by the PB theory because they are highly charged with strong ion correlations.
2.3. Asymmetry problem
It should be noted that there is an inconsistency in the PB equation related with the approximation of $w\left(r\right)\approx q\psi \left(r\right),$ as firstly revealed by Onsager [36]. Let us consider the average concentration of ions of species $i$ at a distance r from a given ion of species $j:$$\begin{eqnarray}{n}_{ji}(r)={n}_{i}\exp (-\beta {w}_{ji}(r)),\end{eqnarray}$and vice versa$\begin{eqnarray}{n}_{ij}(r)={n}_{j}\exp (-\beta {w}_{ij}(r)),\end{eqnarray}$where ${n}_{i}$ is the number density of species $i$ under no constraint. We expect that ${w}_{ij}(r)={w}_{ji}(r)$ since w is the reversible work expended in bringing the two ions $i$ and $j$ from infinity to distance $r.$ The PB theory is on the basis of the approximation that$\begin{eqnarray}{w}_{ji}(r)={q}_{i}{\psi }_{j}(r)\end{eqnarray}$according to equation (14) and$\begin{eqnarray}{w}_{ij}(r)={q}_{j}{\psi }_{i}(r)\end{eqnarray}$according to equation (15). One may expect that$\begin{eqnarray}{q}_{i}{\psi }_{j}(r)={q}_{j}{\psi }_{i}(r).\end{eqnarray}$
The above cannot be true except when $i=j.$ Therefore, the PB solution cannot hold due to nonlinearity except in the symmetric case.
However, it is easy to show that the Yukawa potential in equation (12), the solution for the DH equation, satisfies equation (18). In fact, the DH equation is practically more frequently used because it is linear and has analytic solutions. The most successful application of the DH equation is the calculation of activity coefficient ${\gamma }_{\pm },$ which agrees very well with experiment at the low concentration limit [17],$\begin{eqnarray}\mathrm{ln}\,{\gamma }_{\pm }=-\left|{q}_{+}{q}_{-}\right|\displaystyle \frac{\kappa }{8\pi {\varepsilon }_{0}\varepsilon {k}_{{\rm{B}}}T}.\end{eqnarray}$
Therefore, the DH theory is called the exact limiting law for ionic solutions. Although it is originated from error cancellation, this result suggests that the DH theory is a necessary step in deriving the exact limiting theory. Therefore, the Yukawa potential can be used to develop mean field descriptions for ionic solutions. The DIT [37, 38] developed by Kjellander and Mitchell has an asymptotic solution taking the form of the Yukawa potential but including effective charges and a dielectric constant, which is exact as long as the effective parameters are correctly determined. The MDH theory [39, 40] developed by Xiao and Song suggests a multi-Yukawa potential, which is also exact if the coefficients can be correctly obtained.
3. FT approach
The FT approach provides a way for solving the many-body partition function by transforming the partition function into a functional integration [41–43] through the Hubbard–Stratonovich (HS) transformation [44]. The functional form of the partition function can then be solved by various methods. For example, the simplest approximation is known as the saddle-point approximation which recovers the original PB equation; it can also be solved with the perturbation method [45] and the variational methods [46–50] to obtain more accurate self-consistent equations for depicting ionic systems.
3.1. Formulism
The HS transformation is the key of the FT approach. To understand it, we consider a toy model with the grand canonical partition function:$\begin{eqnarray}{\rm{\Xi }}=\displaystyle \sum _{N=0}^{\infty }\displaystyle \frac{{\lambda }_{{\rm{C}}}^{N}}{N!}{{\rm{e}}}^{-{\varepsilon }_{{\rm{C}}}{N}^{2}/2},\end{eqnarray}$where ${\lambda }_{{\rm{C}}}$ and ${\varepsilon }_{{\rm{C}}}$ are just two mathematical parameters. In order to solve the partition function, one can transform the summation into integrals by noting the two identities$\begin{eqnarray}{{\rm{e}}}^{-{\varepsilon }_{{\rm{C}}}{N}^{2}/2}=\displaystyle \frac{1}{\sqrt{2\pi \varepsilon }}\displaystyle {\int }_{-\infty }^{\infty }{\rm{d}}x{{\rm{e}}}^{-{x}^{2}/2{\varepsilon }_{{\rm{C}}}}{{\rm{e}}}^{{\rm{i}}xN},\end{eqnarray}$$\begin{eqnarray}{{\rm{e}}}^{x}=\displaystyle \sum _{N=0}^{\infty }\displaystyle \frac{{x}^{N}}{N!}.\end{eqnarray}$
The partition function can then be written as$\begin{eqnarray}\begin{array}{l}{\rm{\Xi }}=\displaystyle \sum _{N=0}^{\infty }\displaystyle \frac{{\lambda }_{{\rm{C}}}^{N}}{N!}{{\rm{e}}}^{-{\varepsilon }_{{\rm{C}}}{N}^{2}/2}\\ \,=\displaystyle \frac{1}{\sqrt{2\pi {\varepsilon }_{{\rm{C}}}}}\displaystyle {\int }_{-\infty }^{\infty }{\rm{d}}x{{\rm{e}}}^{-{x}^{2}/2{\varepsilon }_{{\rm{C}}}}\displaystyle \sum _{N=0}^{\infty }\displaystyle \frac{{\lambda }_{{\rm{C}}}^{N}}{N!}{{\rm{e}}}^{{\rm{i}}xN}\\ \,=\displaystyle \frac{1}{\sqrt{2\pi {\varepsilon }_{{\rm{C}}}}}\displaystyle {\int }_{-\infty }^{\infty }{\rm{d}}x\exp \left(-\left({x}^{2}/2{\varepsilon }_{{\rm{C}}}-{\lambda }_{{\rm{C}}}{{\rm{e}}}^{{\rm{i}}x}\right)\right).\end{array}\end{eqnarray}$
The Gaussian identity can be generalized to multiple species as$\begin{eqnarray}{{\rm{e}}}^{-\displaystyle \frac{1}{2}{\varepsilon }_{{\rm{C}}ij}{N}_{i}{N}_{j}}\equiv {{\rm{e}}}^{-\displaystyle \frac{1}{2}({{\boldsymbol{N}}}^{{\bf{T}}}{{\boldsymbol{\varepsilon }}}_{{\rm{C}}}{\boldsymbol{N}})}=\displaystyle {\int }_{-\infty }^{\infty }\displaystyle \frac{\displaystyle {\prod }_{i=1}^{M}{\rm{d}}{x}_{i}}{\sqrt{\det 2\pi {{\boldsymbol{\varepsilon }}}_{{\rm{C}}}}}{{\rm{e}}}^{-\displaystyle \frac{1}{2}({{\boldsymbol{x}}}^{{\bf{T}}}{{{\boldsymbol{\varepsilon }}}_{{\rm{C}}}}^{-1}{\boldsymbol{x}})}{{\rm{e}}}^{{\rm{i}}{{\boldsymbol{N}}}^{{\bf{T}}}{\boldsymbol{x}}},\end{eqnarray}$and to the continuous case [41] as$\begin{eqnarray}\begin{array}{l}{{\rm{e}}}^{-\displaystyle \frac{1}{2}\displaystyle \int {\rm{d}}r{\rm{d}}r^{\prime} \rho \left(r\right){\varepsilon }_{{\rm{C}}}\left(r,r^{\prime} \right)\rho \left(r^{\prime} \right)}\\ \,=\displaystyle \int \displaystyle \frac{{D}\xi }{\sqrt{\det 2\pi {{\boldsymbol{\varepsilon }}}_{{\rm{C}}}}}{{\rm{e}}}^{-\displaystyle \frac{1}{2}\displaystyle \int {\rm{d}}r{\rm{d}}r^{\prime} \xi \left(r\right){{\varepsilon }_{{\rm{C}}}}^{-1}\left(r,r^{\prime} \right)\xi \left(r^{\prime} \right)}{{\rm{e}}}^{{\rm{i}}\displaystyle \int {\rm{d}}r\rho \left(r\right)\xi \left(r\right)},\end{array}\end{eqnarray}$where $\xi \left(r\right)$ is a field variable used to decouple the quadratic interaction on the left hand side.
Let us consider a typical charged particle system consisting of ${n}_{+}$ cations and ${n}_{-}$ anions whose Hamiltonian is$\begin{eqnarray}H=\displaystyle \frac{{e}^{2}}{2}\displaystyle \int {\rm{d}}r{\rm{d}}r^{\prime} \rho \left(r\right)C\left(r,r^{\prime} \right)\rho \left(r^{\prime} \right),\end{eqnarray}$where $e\rho \left(r\right)$ is the total charge density and $C\left(r,r^{\prime} \right)$ is the Coulomb operator defined by$\begin{eqnarray}-{\rm{\nabla }}\cdot \left[{\varepsilon }_{0}\varepsilon {\rm{\nabla }}C\left(r,r^{\prime} \right)\right]=\delta \left(r,r^{\prime} \right).\end{eqnarray}$
The grand partition function of a charged particle system is$\begin{eqnarray}{\rm{\Xi }}=\displaystyle \sum _{{n}_{+}=0}^{\infty }\displaystyle \sum _{{n}_{-}=0}^{\infty }Q\left({n}_{+},{n}_{-}\right){{\rm{e}}}^{{n}_{+}{\mu }_{+}}{{\rm{e}}}^{{n}_{-}{\mu }_{-}},\end{eqnarray}$where ${n}_{\pm }$ are the number densities of cations and anions, ${\mu }_{\pm }$ are the chemical potentials, and $Q\left({n}_{+},{n}_{-}\right)$ is the canonical partition function:$\begin{eqnarray}Q\left({n}_{+},{n}_{-}\right)=\displaystyle \frac{1}{{n}_{+}!{n}_{-}!}\displaystyle \int \displaystyle \prod _{i=1}^{{n}_{+}}{\rm{d}}{r}_{i}\displaystyle \prod _{j=1}^{{n}_{-}}{\rm{d}}{r}_{j}\exp \left(-\beta H\right).\end{eqnarray}$
Applying the HS transformation equation (25) and define a dimensionless field $\varphi \equiv \beta e\xi ,$ the grand partition function becomes [47]$\begin{eqnarray}{\rm{\Xi }}=\displaystyle \frac{1}{{Z}_{{\rm{C}}}}\displaystyle \int {D}\varphi {{\rm{e}}}^{\left\{-L\left[\varphi \right]\right\}},\end{eqnarray}$where ${Z}_{{\rm{C}}}$ is the associated partition function. The action $L\left[\varphi \right]$ in the functional integral is defined as$\begin{eqnarray}L\left[\varphi \right]=\displaystyle \int {\rm{d}}r\left[\displaystyle \frac{1}{2}\in \varphi \left(r\right)\left(-{{\rm{\nabla }}}^{2}\right)\varphi \left(r\right)-{\lambda }_{+}{{\rm{e}}}^{-{\rm{i}}{z}_{+}\varphi }-{\lambda }_{-}{{\rm{e}}}^{-{\rm{i}}{z}_{-}\varphi }\right],\end{eqnarray}$where $\epsilon \equiv {\varepsilon }_{0}\varepsilon /\left(\beta {e}^{2}\right)$ is the rescaled permittivity, ${z}_{\pm }$ are the valences of ions, and ${\lambda }_{\pm }$ are the fugacities [47]. The saddle point for the action, ${\left.\tfrac{\delta L}{\delta \varphi }\right|}_{\varphi ={\varphi }_{0}}=0,$ can be obtained by the Euler–Lagrange equation$\begin{eqnarray}-\epsilon {{\rm{\nabla }}}^{2}{\varphi }_{0}={\rm{i}}{\lambda }_{+}{{\rm{e}}}^{-{\rm{i}}{z}_{+}{\varphi }_{0}}+{\rm{i}}{\lambda }_{-}{{\rm{e}}}^{-{\rm{i}}{z}_{-}{\varphi }_{0}},\end{eqnarray}$which retrogresses into the PB equation if ${\rm{i}}{\varphi }_{0}$ is regarded as the average potential.
3.2. Approximate solutions
The grand partition function equation (30) can be solved with more sophisticated methods to obtain a MPB equation including fluctuation and correlation effects. The variational methods are used by several groups to obtain MPB equations. The most commonly used reference action has the Gaussian form [47, 48, 51]:$\begin{eqnarray}\begin{array}{l}{L}_{{\rm{ref}}}\left[\varphi \right]=\displaystyle \frac{1}{2}\displaystyle \int {\rm{d}}r{\rm{d}}r^{\prime} \left[\varphi \left(r\right)+{\rm{i}}\phi \left(r\right)\right]{G}^{-1}\left(r,r^{\prime} \right)\left[\varphi \left(r^{\prime} \right)\right.\\ \,+\left.{\rm{i}}\phi \left(r^{\prime} \right)\right],\end{array}\end{eqnarray}$where $\varphi \left(r\right)$ is the dimensionless field appeared in equation (30), $G\left(r,r^{\prime} \right)$ is a correlation function which could be the incremental potential in section 4.1, and $\phi \left(r\right)\equiv \beta e\psi \left(r\right)$ corresponds to the average electrostatic field. The variational free energy is$\begin{eqnarray}F={F}_{{\rm{ref}}}+{\left\langle L\left[\varphi \right]-{L}_{{\rm{ref}}}\left[\varphi \right]\right\rangle }_{{\rm{ref}}},\end{eqnarray}$where the average is taken in the reference ensemble with action ${L}_{{\rm{ref}}}.$ Minimizing the free energy functional with respect to $\phi $ and $G,$ we obtain the self-consistent equations [47, 52]$\begin{eqnarray}-{\rm{\nabla }}\cdot \left(\in {\rm{\nabla }}\phi \right)={\lambda }_{+}{z}_{+}{{\rm{e}}}^{-{z}_{+}\phi -{u}_{+}}+{\lambda }_{-}{z}_{-}{{\rm{e}}}^{-{z}_{-}\phi -{u}_{-}},\end{eqnarray}$$\begin{eqnarray}\begin{array}{l}-{\rm{\nabla }}\cdot \left[\in {\rm{\nabla }}G\left(r,r^{\prime} \right)\right]+\left({\lambda }_{+}{z}_{+}^{2}{{\rm{e}}}^{-{z}_{+}\phi -{u}_{+}}+{\lambda }_{-}{z}_{-}^{2}{{\rm{e}}}^{-{z}_{-}\phi -{u}_{-}}\right)\\ \,\times G\left(r,r^{\prime} \right)=\delta \left(r,r^{\prime} \right),\end{array}\end{eqnarray}$$\begin{eqnarray}{u}_{\pm }\left(r\right)=\displaystyle \frac{{z}_{\pm }^{2}}{2}\mathop{\mathrm{lim}}\limits_{r^{\prime} \to r}\left[G\left(r,r^{\prime} \right)-\displaystyle \frac{1}{4\pi \in \left|r-r^{\prime} \right|}\right].\end{eqnarray}$
The term ${u}_{\pm }\left(r\right)$ appeared in the exponent is called the self-energy that includes the correlation and fluctuation effects.
The continuous FT approach is generalized to a hardcore-ion model by Loubet et al [49] to study the liquid–vapor phase transformation. A dipolar PB theory [53] is also developed by treating water molecules as explicit dipoles rather than implicit continuous medium.
The FT approach is useful in calculating dielectric constants [54] but the validity of the Gaussian functional is questionable [51]. We have previously shown that the RDFs calculated by the self-consistent equations, equations (35)–(37), match the molecular dynamics (MD) simulation results only for counter-ions [21]. In addition to the variational methods, a systematic perturbative expansion [41] has been performed by Netz and Orland to solve the partition function by a loop-wise expansion of the action around the saddle point, which indicates that the zero-loop equation is just the ordinary PB equation and the one-loop correction yields similar results as the variational methods.
4. MPB equations with correlations
The ion number distribution in the PB theory is assumed to satisfy the Boltzmann distribution$\begin{eqnarray}\rho \left(r\right)={\rho }_{0}\exp \left(-\beta q\psi \left(r\right)\right),\end{eqnarray}$where ${\rho }_{0}$ is the bulk number density and $q=ze$ is the amount of charge.
Since the Boltzmann distribution is exact only for ideally non-interacting particles, the difference between $q\psi \left(r\right)$ and PMF $w\left(r\right)$ reflects the physical features neglected in the PB theory, i.e. ion size, molecular details of solvent, fluctuation, and ion correlation. On the other hand, the ion distribution with respect to a reference particle, i.e. the RDF, has the following exact relationship with the PMF (equation (8)):$\begin{eqnarray}g\left(r\right)={\rho }_{0}\exp \left(-\beta w\left(r\right)\right).\end{eqnarray}$
Therefore, a more accurate theory may be achieved by dealing with the PMF. In this section, we introduce a correlation-enhanced PB equation developed by us [21], as well as the more sophisticated MPB theories developed by Outhwaite and Bhuiyan [55–57].
4.1. Correlation-enhanced PB model
Considering a two-species ionic system with valences ${z}_{\pm }$ and an average electrostatic potential $\psi \left(r\right),$ the correlation effect can be quantified by inserting a test ion. We initially assume that the system obeys the PB equation, so that when the system is fully relaxed to a new equilibrium state after the insertion of the test ion, the average potential is perturbed to be $\psi \left(r\right)+G\left(r,r^{\prime} \right),$ where $G\left(r,r^{\prime} \right)$ is the incremental potential due to the test ion. For simplicity, we define the dimensionless potential $\phi \equiv \beta e\psi ,$ exactly the same as the one defined in section 3.2, and the Bjerrum length ${l}_{{\rm{B}}}\equiv \tfrac{\beta {e}^{2}}{4\pi {\varepsilon }_{0}\varepsilon },$ then the PB equation with and without the test ion is$\begin{eqnarray}\begin{array}{l}-{{\rm{\nabla }}}^{2}\left[{\phi }_{{\rm{PB}}}(r)+{G}_{\pm }(r,r^{\prime} )\right]=4\pi {l}_{{\rm{B}}}\left[{z}_{+}{\rho }_{{\rm{s}}+}{{\rm{e}}}^{-{z}_{+}{\phi }_{{\rm{PB}}}(r)-{z}_{+}{G}_{\pm }(r,r^{\prime} )}\right.\\ \,+\left.{z}_{-}{\rho }_{{\rm{s}}-}{{\rm{e}}}^{-{z}_{-}{\phi }_{{\rm{PB}}}(r)-{z}_{-}{G}_{\pm }(r,r^{\prime} )}+{z}_{\pm }\delta (r-r^{\prime} )\right]\end{array}\end{eqnarray}$and$\begin{eqnarray}-{{\rm{\nabla }}}^{2}{\phi }_{{\rm{PB}}}=4\pi {l}_{{\rm{B}}}\left({z}_{+}{\rho }_{{\rm{s}}+}{{\rm{e}}}^{-{z}_{+}{\phi }_{{\rm{PB}}}}+{z}_{-}{\rho }_{{\rm{s}}-}{{\rm{e}}}^{-{z}_{-}{\phi }_{{\rm{PB}}}}\right),\end{eqnarray}$where the subscript ‘PB’ stands for the solution to the original PB equation, ${\rho }_{s+}$ and ${\rho }_{s-}$ are bulk number densities of the two ion species. By comparing the two equations, we obtain the expression for $G\left(r,r^{\prime} \right):$$\begin{eqnarray}\begin{array}{l}-{{\rm{\nabla }}}^{2}{G}_{\pm }\left(r,r^{\prime} \right)=4\pi {l}_{{\rm{B}}}\left[{z}_{+}{\rho }_{s+}{{\rm{e}}}^{-{z}_{+}{\phi }_{{\rm{PB}}}\left(r\right)}\left({{\rm{e}}}^{-{z}_{+}{G}_{\pm }\left(r,r^{\prime} \right)}-1\right)\right.\\ \,+\left.{z}_{-}{\rho }_{s-}{{\rm{e}}}^{-{z}_{-}{\phi }_{{\rm{P}}B}\left(r\right)}\left({{\rm{e}}}^{-{z}_{-}{G}_{\pm }\left(r,r^{\prime} \right)}-1\right)+{z}_{\pm }\delta \left(r-r^{\prime} \right)\right].\end{array}\end{eqnarray}$
Equation (42) indicates that two factors contribute to ${G}_{\pm }\left(r,r^{\prime} \right):$ the presence of the test ion and the change of surrounding ions due to the presence of the test ion. Therefore, the correlation effect caused by the insertion of the test ion can be quantified by subtracting $G\left(r,r^{\prime} \right)$ from the potential of the test ion itself ${G}_{0}\left(r,r^{\prime} \right),$ and the correlation term $u\left(r\right),$ which is just the self-energy, is determined by$\begin{eqnarray}{u}_{\pm }(r^{\prime} )=\mathop{\mathrm{lim}}\limits_{r\to r^{\prime} }{u}_{\pm }(r,r^{\prime} )=\mathop{\mathrm{lim}}\limits_{r\to r^{\prime} }\left[{G}_{\pm }(r,r^{\prime} )-{G}_{0\pm }(r,r^{\prime} )\right].\end{eqnarray}$
The self-energy is then included in the original PB equation to obtain the correlation-enhanced PB equation as [21]:$\begin{eqnarray}-{{\rm{\nabla }}}^{2}\phi =4\pi {l}_{{\rm{B}}}\left({z}_{+}{\rho }_{0+}{{\rm{e}}}^{-{z}_{+}\phi -{z}_{+}{u}_{+}}+{z}_{-}{\rho }_{0-}{{\rm{e}}}^{-{z}_{-}\phi -{z}_{-}{u}_{-}}\right),\end{eqnarray}$where ${\rho }_{0\pm }={\rho }_{s\pm }\exp \left({z}_{\pm }{u}_{\infty }\right)$ with ${u}_{\infty }$ being the self-energy at infinity is the rescaled bulk number density to ensure that local density approaches the bulk value when $r$ goes to infinity. As a MPB equation is obtained, the potential ${\phi }_{{\rm{PB}}}$ in equation (42) should be replaced by the solution of equation (44) to obtain the self-consistent equations.
Equations (42)–(44) are similar to the self-consistent equations in the FT approach. It can be shown that the formula of our correlation-enhanced PB model is reduced to the FT approach if one applies a linear approximation to the equation for $G\left(r,r^{\prime} \right)$ [21]. However, such a linearization is in fact mathematically problematic since $G\left(r,r^{\prime} \right)$ diverges as $r^{\prime} $ goes to $r,$ and the linearization is acceptable only when $G\left(r,r^{\prime} \right)$ is close to zero, indicating that the widely used Gaussian reference action is ad hoc and may not be a good choice for charged particles systems. The RDFs obtained by the MD simulations for 1:1 and 1:2 ionic solutions show that our correlation-enhanced PB model matches the MD results very well but the corresponding FT equation deviate significantly from the MD results [21].
4.2. Outhwaite–Bhuiyan MPB theories
Outhwaite and Bhuiyan developed a series of MPB theories known as MPB1 to MPB5 yielding different levels of accuracy [55, 56] with MPB5 the most accurate one. The MPB theories can be considered as a generalization of the integral equation method (see appendix) to the heterogeneous case. The equation describing the model composed of a charged hard-sphere particle next to a planar hard electrode with a uniform surface charge density $\sigma $ is [58]$\begin{eqnarray}\displaystyle \frac{{{\rm{d}}}^{2}\psi }{{\rm{d}}{x}^{2}}=-\displaystyle \frac{1}{{\varepsilon }_{0}\varepsilon }\displaystyle \displaystyle \sum _{s}{q}_{s}{\rho }_{s}{g}_{s}\left(x\right),\end{eqnarray}$where ${g}_{s}\left(x\right)$ is the singlet wall-ion distribution function which can be obtained using the Kirkwood charging process [59, 60] by$\begin{eqnarray}{g}_{s}\left(x\right)={\xi }_{s}\left(x\right)\exp \left[-\beta {q}_{s}\psi \left(x\right)-\beta \displaystyle {\int }_{0}^{{e}_{s}}\mathop{\mathrm{lim}}\limits_{r\to 0}\left({\varphi }_{s}^{* }-{\varphi }_{sb}^{* }\right){\rm{d}}{q}_{s}\right],\end{eqnarray}$where ${\xi }_{s}\left(x\right)={g}_{s}\left(x\left|{q}_{s}=0\right.\right)$ is the exclusive volume term, ${\varphi }^{* }$ is the fluctuation potential, the subscript $b$ denotes the corresponding value for bulk, and the subscript $s$ denotes the ion species. After complicated calculation of the fluctuation potential, the expression finally appears to be$\begin{eqnarray}{g}_{s}\left(x\right)={\xi }_{s}\left(x\right)\exp \left[-\beta {q}_{s}L\left(\psi \right)-\displaystyle \frac{\beta {q}^{2}}{8\pi {\varepsilon }_{0}\varepsilon }\left(F-{F}_{0}\right)\right],\end{eqnarray}$where the operator $L\left(\psi \right)$ is given by$\begin{eqnarray}\begin{array}{l}L\left(\psi \right)=\displaystyle \frac{F}{2}\left\{\psi \left(x+a\right)+\psi \left(x-a\right)\right\}\\ \,-\displaystyle \frac{\left(F-{F}_{0}\right)}{2a}\displaystyle {\int }_{x-a}^{x+a}\psi \left(X\right){\rm{d}}X,\end{array}\end{eqnarray}$and the coefficients $F$ and ${F}_{0}$ are expressed in terms of the distance $x$ and ion diameter $a.$ Details of the derivation can be found in [56]. The expression for the exclusive volume term can be obtained by applying the Born–Green–Yvon hierarchy [61], or by using an alternative and simpler expression developed by Outhwaite and Lamperski [62]:$\begin{eqnarray}{\xi }_{s}\left(x\right)={\xi }_{s}^{0}\left(x\right)\exp \left[\displaystyle \int \displaystyle \displaystyle \sum _{t}{\rho }_{t}{c}_{st}^{0}\left\{{g}_{t}\left(x\right)-{g}_{t}\left(x\left|\sigma =0\right.\right)\right\}{\rm{d}}V\right],\end{eqnarray}$where ${c}_{st}^{0}$ is the inhomogeneous neutral-sphere direct correlation function (DCF) whose meaning can be found in the appendix and ${\xi }_{s}^{0}\left(x\right)={\xi }_{s}\left(x\left|\sigma =0\right.\right).$ The ${c}_{st}^{0}$ and ${\xi }_{s}^{0}\left(x\right)$ can be approximated by the PY results [63, 64] and the subscript $t$ runs over all the ion species.
The MPB equations are solved for a primitive model of a planar double layer containing 1:1 or 1:2 charged hard-sphere particles. The results are compared with corresponding Monte Carlo (MC) simulation results and both the MPB formulations reproduce the MC data very well [58].
5. Mean field theories
The MPB equations described above are usually hard or impossible to be solved analytically, which limits their practical applications in some scenarios. A simple analytic expression for the ion or potential distribution is important and convenient in many applications. In this section, we introduce two mean-field theories, the DIT [37, 38, 65] and the MDH theory [39, 40]. The former provides an effective-charge Yukawa potential and the latter ends up with a multi-Yukawa potential.
5.1. The DIT
To derive the DIT, we start with the Ornstein–Zernike (OZ) equation for the pair correlation functions:$\begin{eqnarray}{h}_{ij}\left(r\right)={c}_{ij}\left(r\right)+\displaystyle \displaystyle \sum _{l}\displaystyle \int {\rm{d}}r^{\prime} {c}_{il}\left(\left|r-r^{\prime} \right|\right){n}_{l}{h}_{lj}\left(r^{\prime} \right),\end{eqnarray}$where $h\left(r\right)$ is the total correlation function (TCF) and $c\left(r\right)$ is the DCF, see the appendix for their descriptions. It is convenient to apply the Fourier transform to obtain$\begin{eqnarray}{\hat{h}}_{ij}\left(k\right)={\hat{c}}_{ij}\left(k\right)+\displaystyle \displaystyle \sum _{l}{\hat{c}}_{il}\left(k\right){n}_{l}{\hat{h}}_{lj}\left(k\right)\end{eqnarray}$and convert it to the matrix form$\begin{eqnarray}\hat{{\boldsymbol{H}}}=\hat{{\boldsymbol{C}}}+\hat{{\boldsymbol{C}}}N\hat{{\boldsymbol{H}}},\end{eqnarray}$where $\hat{{\boldsymbol{H}}}=\left\{{\hat{h}}_{ij}\right\},\hat{{\boldsymbol{C}}}=\left\{{\hat{c}}_{ij}\right\},\,{\boldsymbol{N}}=\left\{{\delta }_{ij}{n}_{i}\right\}.$ The total charge distribution is$\begin{eqnarray}{\rho }_{j}^{{\rm{tot}}}\left(r\right)={q}_{j}{\delta }^{\left(3\right)}\left(r\right)+\displaystyle \displaystyle \sum _{l}{q}_{l}{n}_{l}{h}_{lj}\left(r\right)\end{eqnarray}$and the average potential is$\begin{eqnarray}{\psi }_{j}\left(r\right)=\displaystyle \int {\rm{d}}r^{\prime} \varphi \left(\left|r-r^{\prime} \right|\right){\rho }_{j}^{{\rm{tot}}}\left(r^{\prime} \right),\end{eqnarray}$where $\varphi \left(r\right)=1/\left(4\pi \varepsilon {\varepsilon }_{0}r\right)$ is the Coulomb potential. In the Fourier space the average potential is$\begin{eqnarray}{\hat{\psi }}_{j}\left(k\right)=\hat{\varphi }\left(k\right){\hat{\rho }}_{j}^{{\rm{tot}}}\left(k\right)=\hat{\varphi }\left(k\right)\left[{q}_{j}+\displaystyle \displaystyle \sum _{l}{q}_{l}{n}_{l}{\hat{h}}_{lj}\left(k\right)\right]\end{eqnarray}$and in the matrix form$\begin{eqnarray}{\left(\hat{{\boldsymbol{\psi }}}\right)}^{{\rm{T}}}=\hat{\varphi }\left[{{\boldsymbol{q}}}^{{\rm{T}}}+{{\boldsymbol{q}}}^{{\rm{T}}}N\hat{{\boldsymbol{H}}}\right],\end{eqnarray}$where $\hat{{\boldsymbol{\psi }}}=\left\{{\hat{\psi }}_{j}\right\},{\boldsymbol{q}}=\left\{{q}_{j}\right\}$ are vectors.
The key idea of the DIT is to split the pair correlation function into a short-range part and a long-range part, then the DCF is$\begin{eqnarray}{c}_{ij}\left(r\right)={c}_{ij}^{0}\left(r\right)+{c}_{ij}^{l}\left(r\right),\end{eqnarray}$where superscript 0 denotes the short-range part and superscript l denotes the long-range part. Analogous to the long-range behavior of the RDF, the long-range part of DCF ${c}_{ij}^{l}\left(r\right)\sim -\beta {q}_{i}{q}_{j}\varphi \left(r\right)$ when $r\to \infty .$ Therefore, we have$\begin{eqnarray}{C}_{ij}^{0}\left(r\right)={C}_{ij}\left(r\right)+\beta {q}_{i}{q}_{j}\varphi \left(r\right)\end{eqnarray}$and the Fourier transform in the matrix form is$\begin{eqnarray}{\hat{{\boldsymbol{C}}}}^{0}=\hat{{\boldsymbol{C}}}+\beta {\boldsymbol{q}}{{\boldsymbol{q}}}^{{\rm{T}}}\hat{\varphi }.\end{eqnarray}$
Inserting the short-range part of DCF into the OZ equation and utilizing equation (56), we have$\begin{eqnarray}\hat{{\boldsymbol{H}}}={\hat{{\boldsymbol{C}}}}^{0}+{\hat{{\boldsymbol{C}}}}^{0}N\hat{{\boldsymbol{H}}}-\beta {\boldsymbol{q}}{\left(\hat{{\boldsymbol{\psi }}}\right)}^{{\rm{T}}}.\end{eqnarray}$
The solution of the TCF is$\begin{eqnarray}\hat{{\boldsymbol{H}}}={\left({\bf{1}}-{\hat{{\boldsymbol{C}}}}^{0}{\boldsymbol{N}}\right)}^{-1}{\hat{{\boldsymbol{C}}}}^{0}-\beta {\left({\bf{1}}-{\hat{{\boldsymbol{C}}}}^{0}{\boldsymbol{N}}\right)}^{-1}{\boldsymbol{q}}{\left(\hat{{\boldsymbol{\psi }}}\right)}^{{\rm{T}}}.\end{eqnarray}$
Defining the short-range variables ${\hat{{\boldsymbol{H}}}}^{0}={\left({\bf{1}}-{\hat{{\boldsymbol{C}}}}^{0}{\boldsymbol{N}}\right)}^{-1}{\hat{{\boldsymbol{C}}}}^{0}$ and ${\hat{\rho }}^{0}={\left({\bf{1}}-{\hat{{\boldsymbol{C}}}}^{0}{\boldsymbol{N}}\right)}^{-1}{\boldsymbol{q}},$ we have$\begin{eqnarray}\hat{{\boldsymbol{H}}}={\hat{{\boldsymbol{H}}}}^{0}-\beta {\hat{{\boldsymbol{\rho }}}}^{0}{\left(\hat{{\boldsymbol{\psi }}}\right)}^{{\rm{T}}}.\end{eqnarray}$
In the real space, the TCF is$\begin{eqnarray}{h}_{ij}\left(r\right)={h}_{ij}^{0}\left(r\right)-\beta \displaystyle \int {\rm{d}}r^{\prime} {\rho }_{i}^{0}\left(\left|r-r^{\prime} \right|\right){\psi }_{j}\left(r^{\prime} \right).\end{eqnarray}$
The form of equation (63) is the same as that of the DH theory, but their definitions of the charge densities are different. In the DH theory, the corresponding expression for the TCF is$\begin{eqnarray}{h}_{ij}^{{\rm{DH}}}\left(r\right)=-\beta \displaystyle \int {\rm{d}}r^{\prime} {q}_{i}{\delta }^{\left(3\right)}\left(r\right){\psi }_{j}^{{\rm{DH}}}\left(r^{\prime} \right),\end{eqnarray}$where ${\delta }^{\left(3\right)}\left(r\right)$ is the three-dimensional Dirac delta function. For DIT, to derive an expression for the average potential, we first express ${\rho }^{0}$ in terms of ${h}_{ij}:$$\begin{eqnarray}{\hat{{\boldsymbol{\rho }}}}^{0}={\boldsymbol{q}}+{\hat{{\boldsymbol{H}}}}^{0}{\boldsymbol{Nq}}.\end{eqnarray}$
In the real space$\begin{eqnarray}{\rho }_{j}^{0}\left(r\right)={q}_{j}{\delta }^{\left(3\right)}\left(r\right)+\displaystyle \displaystyle \sum _{l}{q}_{l}{n}_{l}{h}_{lj}^{0}\left(r\right).\end{eqnarray}$
Combining equations (53), (63), and (66), we obtain the total density$\begin{eqnarray}\begin{array}{l}{\rho }_{j}^{{\rm{tot}}}\left(r\right)={q}_{j}{\delta }^{\left(3\right)}\left(r\right)+\displaystyle \displaystyle \sum _{l}{q}_{l}{n}_{l}{h}_{lj}\left(r\right)\\ \,={\rho }_{j}^{0}\left(r\right)-\beta \displaystyle \int {\rm{d}}r^{\prime} \displaystyle \displaystyle \sum _{l}{q}_{l}{n}_{l}{\rho }_{l}^{0}\left(\left|r-r^{\prime} \right|\right){\psi }_{j}\left(r^{\prime} \right)\\ \,={\rho }_{j}^{0}\left(r\right)-\displaystyle \int {\rm{d}}r^{\prime} \alpha \left(\left|r-r^{\prime} \right|\right){\psi }_{j}\left(r^{\prime} \right),\end{array}\end{eqnarray}$where $\alpha \left(\left|r-r^{\prime} \right|\right)\equiv \displaystyle {\sum }_{l}{q}_{l}{n}_{l}{\rho }_{l}^{0}\left(\left|r-r^{\prime} \right|\right)$ is a response function comparable to the Debye parameter ${\kappa }_{{\rm{DH}}}^{2}=\beta \displaystyle {\sum }_{l}{n}_{l}{q}_{l}^{2}.$
The average potential satisfies Poisson’s equation:$\begin{eqnarray}-{\varepsilon }_{0}\varepsilon {{\rm{\nabla }}}^{2}{\psi }_{j}\left(r\right)={\rho }_{j}^{{\rm{tot}}}\left(r\right)={\rho }_{j}^{0}\left(r\right)-\displaystyle \int {\rm{d}}r^{\prime} \alpha \left(\left|r-r^{\prime} \right|\right){\psi }_{j}\left(r^{\prime} \right),\end{eqnarray}$which reduces to the DH equation if we replace ${\rho }_{j}^{0}\left(r\right)$ with the ion density and $\alpha \left(\left|r-r^{\prime} \right|\right)$ with ${\kappa }_{{\rm{DH}}}^{2}.$ We see that all information we need to obtain the solution is the expression of $\alpha \left(\left|r-r^{\prime} \right|\right).$ Applying the linear response theory [66], it can be shown that $\alpha \left(r\right)$ is related to the dielectric function ${\hat{\varepsilon }}_{l}\left(k\right)$ in the Fourier space by$\begin{eqnarray}\hat{\varepsilon }\left(k\right)=\varepsilon +\displaystyle \frac{\hat{\alpha }\left(k\right)}{{\varepsilon }_{0}{k}^{2}},\end{eqnarray}$where $\hat{\alpha }\left(k\right)$ is the Fourier transform of $\alpha \left(r\right).$ Therefore, all the information of the key parameters in DIT is contained in the response functions. Practically, the linear response function $\alpha \left(r\right)$ can be obtained by experiment and/or MD simulation. In the theoretical work by Varela et al [67], $\alpha \left(r\right)$ is also obtained by calculating the static structure factor for a one-component charged-sphere model in the modified mean-spherical approximation (MSA) and the random phase approximation (RPA). The authors have shown that the RPA clearly accounts for the basic feature of the DIT that the distribution functions split into two well-defined parts and can be considered as a first-order approximation to the problem [12].
5.2. The MDH theory
The MDH theory, which is closely related to the DIT, still needs the information of the dielectric function to determine its key factors. An advantage of the MDH theory is that the coefficients are not too sensitive to the size and charge of the solute molecules [39], which allows the usage of the solvent to estimate the coefficients for the solute. The Gaussian unit is adopted and the $\unicode{x002C6}$ on top of the variables in the Fourier space is omitted in the following equations, as is in the original paper.
To derive the MDH theory, we start with the Poisson’s equation in the Fourier space ([40]) of the DIT:$\begin{eqnarray}{k}^{2}\varepsilon \left(k\right){\psi }_{j}\left(k\right)=4\pi {\rho }_{j}^{0}\left(k\right),\end{eqnarray}$in which $\varepsilon \left(k\right)=\varepsilon +\tfrac{4\pi \alpha \left(k\right)}{{k}^{2}}$ is the dielectric function of the solution with the relative dielectric constant $\varepsilon $ and $\alpha \left(k\right)$ is just $\hat{\alpha }\left(k\right)$ in equation (69). The charge density ${\rho }_{j}^{0}\left(k\right)$ has the same definition as in equation (66), but it can be interpreted as the free charge density in most cases as long as the linear response theory is applicable [40].
The roots of $\varepsilon \left(k\right)$ are useful for finding the solution and are denoted as ${\rm{i}}{k}_{n},$ therefore$\begin{eqnarray}{k}^{2}\varepsilon \left(k\right)=\displaystyle \frac{\varepsilon \displaystyle {\prod }_{n}\left({k}^{2}+{k}_{n}^{2}\right)}{f\left(k\right)},\end{eqnarray}$where $f\left(k\right)$ is a function with no poles. We can then solve the potential in the real space using the inverse Fourier transform$\begin{eqnarray}\begin{array}{l}{\psi }_{j}\left(r\right)=\displaystyle \frac{1}{2{\pi }^{2}r}\displaystyle {\int }_{0}^{\infty }{\rm{d}}k\left[{\psi }_{j}\left(k\right)k\right]\sin \left(kr\right)\\ \,=\displaystyle \frac{1}{2{\pi }^{2}r}\displaystyle {\int }_{0}^{\infty }{\rm{d}}k\left[\displaystyle \frac{kf\left(k\right){\rho }_{j}^{0}\left(k\right)}{\varepsilon \displaystyle {\prod }_{n}\left({k}^{2}+{k}_{n}^{2}\right)}\right]\sin \left(kr\right).\end{array}\end{eqnarray}$
Solving the integral using the residue theorem, we have the following asymptotic solution:$\begin{eqnarray}{\psi }_{j}\left(r\right)\sim \displaystyle \displaystyle \sum _{n}\displaystyle \frac{{q}_{{\rm{eff}},jn}}{{\varepsilon }_{{\rm{eff}},n}}\displaystyle \frac{{{\rm{e}}}^{-{k}_{n}r}}{r},\end{eqnarray}$where ${q}_{{\rm{eff}},jn}={\rho }_{j}^{0}\left({\rm{i}}{k}_{n}\right)$ is an effective charge and ${\varepsilon }_{{\rm{eff}},n}=\displaystyle \frac{1}{2}{\left[k\tfrac{{\rm{d}}\varepsilon \left(k\right)}{{\rm{d}}k}\right]}_{k={\rm{i}}{k}_{n}}$ is an effective dielectric constant. The solution is a linear combination of multiple Yukawa potentials with different Debye lengths ${k}_{n}.$ The DH theory can be obtained by applying $\varepsilon \left(k\right)=\varepsilon \left(1+{k}_{{\rm{DH}}}^{2}/{k}^{2}\right)$ to the MDH theory so that it has only one root ${\rm{i}}{k}_{{\rm{DH}}}$ and hence results in a Yukawa potential.
Since the solution of the MDH equation is a combination of the solutions of DH equations with various Debye lengths, it can be determined by finding an effective linearly-combined coefficient set instead of using the effective parameters directly. Given the dielectric function is known, this kind of coefficient set can be determined by the Stillinger–Lovett relation [68–70]. In practice, the dielectric function can be calculated by the MSA or the HNC approximation.
6. Conclusions
In this paper, we review the original PB theory and some recent progresses of improving it. The original PB theory is the first successful theoretical model for describing charged systems including ionic solutions, but the neglect of ion correlations and exclusive volume effects disables its applicability to the cases of high charge densities or valences.
The FT approach, aided by the definition of self-energy, provides a way to directly solve the grand partition function by performing the HS transformation, and various approximation methods lead to different versions of the FT approach. The variational method with a Gaussian reference action was used by several groups to solve the functional form of the partition function. However, the validity of the Gaussian functional is questionable, and we found its numerical solution deviates significantly from MD simulation results.
A correlation-enhanced PB model was developed by us based on physical pictures, which also defines a self-energy but whose expression is different from the FT approach. We have shown that our model reduces to the variational method of the FT approach when applying a linear approximation, which is nevertheless mathematically problematic. A more sophisticated MPB theory was developed by Outhwaite and Bhuiyan, which still keeps the form of the original PB equation, but the pair correlation functions are calculated analytically with the integral equation theory.
The mean field theories try to find a solution as simple as the Yukawa potential. The DIT gives an exact expression, which has the form of the DH theory with the variables to be determined by splitting the pair correlation function into a short-range part and a long-range part. The MDH theory suggests a multi-Yukawa potential which is also exact as long as the coefficients are precisely determined. However, both the DIT and MDH theories require the information of dielectric function, which needs to be further determined either by theoretical approximations or experimental data.
Acknowledgments
This work was supported by the Strategic Priority Research Program of Chinese Academy of Sciences (Grant No. XDA17010504), the National Natural Science Foundation of China (Nos. 11774357 and 11947302), and the CAS Biophysics Interdisciplinary Innovation Team Project (No. 2060299).
Appendix. Integral equation theory
The DCF $c\left(r\right)$ and TCF $h\left(r\right)$ are related by the OZ equation:$\begin{eqnarray}{h}_{ij}\left(r\right)={c}_{ij}\left(r\right)+\displaystyle \displaystyle \sum _{l}\displaystyle \int {\rm{d}}r^{\prime} {c}_{il}\left(\left|r-r^{\prime} \right|\right){n}_{l}{h}_{lj}\left(r^{\prime} \right).\end{eqnarray}$
The OZ equation can be understood as follows: the particles $i$ and $j$ can be related directly by ${c}_{ij}\left(r\right),$ or indirectly through another particle $l.$ The total correlation function is defined as$\begin{eqnarray}h\left(r\right)=g\left(r\right)-1,\end{eqnarray}$where $g\left(r\right)$ is the RDF. Once we obtain the total correlation function or equivalently the RDF, the properties of the system can be calculated with the knowledge of statistical physics. However, the DCF is hard to determine, not to mention that its physical meaning is not well defined. To solve the OZ equation, several approximate methods, known as the integral equation theories [24], are proposed.
An intuitive method is the MSA, which is defined in terms of the DCF by:$\begin{eqnarray}\begin{array}{ll}g\left(r\right)=0, & r\lt d\\ c\left(r\right)=-\beta v\left(r\right), & r\gt d\end{array}\end{eqnarray}$where $d$ is the hard sphere diameter and $v\left(r\right)$ is the interaction potential between particles. In the MSA, the OZ equation can be solved analytically and can be shown to reduce to the DH theory at a high temperature or a low density limit.
Other classic approximations are the PY and HNC equations:$\begin{eqnarray}\begin{array}{ll}g\left(r\right)=\exp \left[-\beta v\left(r\right)\right]\left[1+h\left(r\right)-c\left(r\right)\right] & \left({\rm{PY}}\right),\\ g\left(r\right)=\exp \left[-\beta v\left(r\right)\right]\exp \left[h\left(r\right)-c\left(r\right)\right] & \left({\rm{HNC}}\right).\end{array}\end{eqnarray}$
The PY or HNC equation can be solved numerically together with the OZ equation. The HNC approximation can reproduce simulation results for charged particles systems very accurately, while the PY approximation is proved to be more successful for short-range potentials.
KjellanderR2001Distribution Function Theory of Electrolytes and Electrical Double Layers—Charge Renormalisation and Dressed Ion TheoryElectrostatic Effects in Soft Matter and Biophysics vol 46DordrechtSpringer [Cited within: 2]
近日，中国科学院遗传与发育生物学研究所李家洋研究组受邀在权威学术期刊National Science Review发表题为《Short and Long Term Challenges in Crop Breeding》的社论文章(DOI:10.1093/nsr/nwab002)，文中对作物育种所面 ...
,1,21CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, 55 East Zhongguancun Road, PO Box 2735, Beijing 100190, 
New window|Download| PPT slide
