删除或更新信息,请邮件至freekaoyan#163.com(#换成@)

A brief review of continuous models for ionic solutions: the Poisson【-逻*辑*与-】ndash;Boltzmann and rel

本站小编 Free考研考试/2022-01-02

闂傚倸鍊搁崐椋庣矆娴h櫣绀婂┑鐘插€寸紓姘辨喐閺冨牄鈧線寮介鐐茶€垮┑锛勫仧缁垶寮悩缁樷拺闂侇偆鍋涢懟顖涙櫠閹绢喗鐓熼柟鍨暙娴滄壆鈧娲栨晶搴ㄥ箲閸曨剚濮滈柡澶嬪閻庢娊姊婚崒娆戠獢闁逞屽墰閸嬫盯鎳熼娑欐珷濞寸厧鐡ㄩ悡鏇㈡煟濡崵鍙€闁告瑥瀚埀顒冾潐濞插繘宕归懞銉ょ箚闁割偅娲栭悙濠囨煏婵炲灝鍔村ù鍏兼礋濮婃椽鎳¢妶鍛€鹃梺鑽ゅ枂閸庢娊鍩€椤掍礁鍤柛娆忓暙椤曪綁骞庨挊澶愬敹闂侀潧顧€婵″洭宕㈤柆宥嗏拺鐟滅増甯掓禍浼存煕閹惧鎳囬柕鍡楀暙閳诲酣骞嬮悩纰夌床闂備礁鎲¢悷锕傛晪閻庤鎸稿Λ娑㈠焵椤掑喚娼愭繛鎻掔箻瀹曟繈骞嬮敃鈧弸渚€鏌熼崜褏甯涢柡鍛倐閺屻劑鎮ら崒娑橆伓闂傚倸鍊搁崐鐑芥倿閿旈敮鍋撶粭娑樻噽閻瑩鏌熸潏楣冩闁稿孩鏌ㄩ埞鎴﹀磼濮橆厼鏆堥梺绋款儑閸犳劗鎹㈠☉銏犵婵炲棗绻掓禒鑲╃磼缂併垹骞愰柛瀣崌濮婅櫣鎷犻弻銉偓妤佺節閳ь剚娼忛妸锕€寮块柣搴ㄦ涧閹芥粍绋夊鍡愪簻闁哄稁鍋勬禒锕傛煟閹惧瓨绀冪紒缁樼洴瀹曞崬螣閸濆嫷娼曞┑鐘媰鐏炶棄顫紓浣虹帛缁诲牓宕洪埀顒併亜閹烘垵顏╃紒鐘劜閵囧嫰寮埀顒勫磿閸愯尙鏆﹂柕澶堝劗閺€浠嬫煟閹邦剙绾фい銉у仱閹粙顢涘⿰鍐ф婵犵鈧磭鍩fい銏℃礋閺佸倿鎮剧仦钘夌闂傚倷鑳舵灙闁哄牜鍓涚划娆撳箻鐠囪尙鐤囬梺绯曞墲閻燂箓宕戦弽銊х闁糕剝蓱鐏忎即鏌i幘瀛樼濞e洤锕、娑樷枎閹烘繂濡抽梻浣呵圭€涒晠宕归崷顓燁潟闁规崘顕х壕鍏兼叏濡搫鎮戝Δ鏃堟⒒娓氣偓閳ь剛鍋涢懟顖涙櫠鐎涙ḿ绠惧ù锝呭暱閸氭ê鈽夊Ο閿嬵潔闂侀潧绻嗛埀顒€鍘栭崙鑺ョ節閻㈤潧孝闁挎洏鍊濋幃褎绻濋崶褏鏌у銈嗗笒鐎氼參鎮¢妷鈺傜厽闁哄洨鍋涢埀顒€婀遍埀顒佺啲閹凤拷
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

Received:2020-01-22Revised:2020-03-28Accepted:2020-03-29Online:2020-05-22


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 [13] and the structure and dynamics of biomolecules [46]. 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 [1214]. 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)$.

Combining equations (5), (6), (8) together, we have$\begin{eqnarray}{{\rm{\nabla }}}^{2}\left\langle \psi (r)\right\rangle =-\displaystyle \frac{1}{{\varepsilon }_{0}\varepsilon }\displaystyle \displaystyle \sum _{s}{c}_{s}{q}_{s}{{\rm{e}}}^{-\beta {w}_{1s}(r)}.\end{eqnarray}$

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 [2529], like-charge attraction [3034], 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 [46, 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 [4143] 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 [4650] 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 [5557].

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 [6870]. 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.

Reference By original order
By published year
By cited within times
By Impact factor

Wang R Wang Z G 2013 J. Chem. Phys. 139124702
DOI:10.1063/1.4821636 [Cited within: 1]

Wang R Wang Z G 2014 Phys. Rev. Lett. 112136101
DOI:10.1103/PhysRevLett.112.136101

Adar R M Andelman D Diamant H 2016 Phys. Rev. E 94022803
DOI:10.1103/PhysRevE.94.022803 [Cited within: 1]

Tan Z Zhang W Shi Y Wang F 2015 Advance in Structural BioinformaticsWei D et al.DordrechtSpringer143
[Cited within: 2]

Tan Z-J Chen S-J 2005 J. Chem. Phys. 122044903
DOI:10.1063/1.1842059

Tan Z J Chen S J 2006 Biophys. J. 90 1175
DOI:10.1529/biophysj.105.070904 [Cited within: 3]

Debye P Huckel E 1923 Phys. Z. 24 185
[Cited within: 2]

Andelman D 1995 Handbook of Biological PhysicsAmsterdamElsevier603
[Cited within: 1]

Mester Z Panagiotopoulos A Z 2015 J. Chem. Phys. 142044507
DOI:10.1063/1.4906320 [Cited within: 1]

Levin Y 2002 Rep. Prog. Phys. 65 1577
DOI:10.1088/0034-4885/65/11/201 [Cited within: 1]

Honig B Nicholls A 1995 Science 268 1144
DOI:10.1126/science.7761829 [Cited within: 1]

Varela L M Garcia M Mosquera V 2003 Phys. Rep. 382 1
DOI:10.1016/S0370-1573(03)00210-2 [Cited within: 2]

Grochowski P Trylska J 2008 Biopolymers 89 93
DOI:10.1002/bip.20877

Messina R 2009 J. Phys.: Condens. Matter 21113102
DOI:10.1088/0953-8984/21/11/113102 [Cited within: 1]

Li H Lu B 2014 J. Chem. Phys. 141024115
DOI:10.1063/1.4887342 [Cited within: 1]

Liu X Qiao Y Lu B 2018 SIAM J. Appl. Math. 78 1131
DOI:10.1137/16M1108583 [Cited within: 1]

McQuarrie D A 1976 Statistical MechanicsNew YorkHarper and Row
[Cited within: 5]

Chandler D 1987 Introduction to Modern Statistical MechanicsOxfordOxford University Press
[Cited within: 1]

Bolt G H 1955 J. Colloid Sci. 10 206
DOI:10.1016/0095-8522(55)90027-1 [Cited within: 1]

Lu B Z Zhou Y C Holst M J McCammon J A 2008 Commun. Comput. Phys. 3 973
[Cited within: 1]

Su M Xu Z Wang Y 2019 J. Phys.: Condens. Matter 31355101
DOI:10.1088/1361-648X/ab24a9 [Cited within: 6]

Ren G Wang Y 2014 Eur. Phys. Lett. 10730005
DOI:10.1209/0295-5075/107/30005 [Cited within: 1]

Roux B. t. Simonson T 1999 Biophys. Chem. 78 1
DOI:10.1016/S0301-4622(98)00226-9 [Cited within: 1]

Hansen J-P McDonald I R 2006 Theory of Simple Liquids3rd ednHansen J-P McDonald I RNew YorkAcademic78
[Cited within: 2]

Grosberg A Y Nguyen T T Shklovskii B I 2002 Rev. Mod. Phys. 74 329
DOI:10.1103/RevModPhys.74.329 [Cited within: 2]

Shklovskii B I 1999 Phys. Rev. E 60 5802
DOI:10.1103/PhysRevE.60.5802

Gonzalez-Tovar E Jimenez-Angeles F Messina R Lozada-Cassou M 2004 J. Chem. Phys. 120 9782
DOI:10.1063/1.1710861

Messina R Holm C Kremer K 2001 Phys. Rev. E 64021405
DOI:10.1103/PhysRevE.64.021405

Yu J Aguilar-Pineda G E Antillon A Dong S H Lozada-Cassou M 2006 J. Colloid Interface Sci. 295 124
DOI:10.1016/j.jcis.2005.08.016 [Cited within: 1]

Diehl A Carmona H A Levin Y 2001 Phys. Rev. E 64011804
DOI:10.1103/PhysRevE.64.011804 [Cited within: 1]

Gronbech-Jensen N Mashl R J Bruinsma R F Gelbart W M 1997 Phys. Rev. Lett. 78 2477
DOI:10.1103/PhysRevLett.78.2477

Ha B Y Liu A J 1997 Phys. Rev. Lett. 79 1289
DOI:10.1103/PhysRevLett.79.1289

Kardar M Golestanian R 1999 Rev. Mod. Phys. 71 1233
DOI:10.1103/RevModPhys.71.1233

Podgornik R Parsegian V A 1998 Phys. Rev. Lett. 80 1560
DOI:10.1103/PhysRevLett.80.1560 [Cited within: 1]

Lin C Zhang X Qiang X Zhang J-S Tan Z-J 2019 J. Chem. Phys. 151114902
DOI:10.1063/1.5120756 [Cited within: 1]

Onsager L 1933 Chem. Rev. 13 73
DOI:10.1021/cr60044a006 [Cited within: 1]

Kjellander R 2001 Distribution 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]

Kjellander R Mitchell D J 1994 J. Chem. Phys. 101 603
DOI:10.1063/1.468116 [Cited within: 2]

Xiao T 2015 ChemPhysChem 16 833
DOI:10.1002/cphc.201402694 [Cited within: 3]

Xiao T Song X 2011 J. Chem. Phys. 135104104
DOI:10.1063/1.3632052 [Cited within: 4]

Frydel D 2015 Eur. J. Phys. 36065050
DOI:10.1088/0143-0807/36/6/065050 [Cited within: 3]

Netz R R Orland H 1999 Europhys. Lett. 45 726
DOI:10.1209/epl/i1999-00228-6

Naji A Kanduc M Forsman J Podgornik R 2013 J. Chem. Phys. 139 13
DOI:10.1063/1.4824681 [Cited within: 1]

Hubbard J 1959 Phys. Rev. Lett. 3 77
DOI:10.1103/PhysRevLett.3.77 [Cited within: 1]

Netz R R Orland H 2000 Eur. Phys. J. E 1 203
DOI:10.1007/s101890050023 [Cited within: 1]

Netz R R Orland H 2003 Eur. Phys. J. E 11 301
DOI:10.1140/epje/i2002-10159-0 [Cited within: 1]

Wang Z G 2010 Phys. Rev. E 81021501
DOI:10.1103/PhysRevE.81.021501 [Cited within: 4]

Xu Z L Maggs A C 2014 J. Comput. Phys. 275 310
DOI:10.1016/j.jcp.2014.07.004 [Cited within: 1]

Loubet B Manghi M Palmeri J 2016 J. Chem. Phys. 145044107
DOI:10.1063/1.4959034 [Cited within: 1]

Buyukdagli S Manghi M Palmeri J 2010 Phys. Rev. E 81041601
DOI:10.1103/PhysRevE.81.041601 [Cited within: 1]

Buyukdagli S Blossey R 2016 J. Phys.: Condens. Matter 28343001
DOI:10.1088/0953-8984/28/34/343001 [Cited within: 2]

Wang R Wang Z G 2015 J. Chem. Phys. 142104705
DOI:10.1063/1.4914170 [Cited within: 1]

Abrashkin A Andelman D Orland H 2007 Phys. Rev. Lett. 99077801
DOI:10.1103/PhysRevLett.99.077801 [Cited within: 1]

Levy A Andelman D Orland H 2012 Phys. Rev. Lett. 108227801
DOI:10.1103/PhysRevLett.108.227801 [Cited within: 1]

Outhwaite C W Bhuiyan L B Levine S 1980 J. Chem. Soc. Faraday Trans. II 76 1388
DOI:10.1039/F29807601388 [Cited within: 2]

Outhwaite C W Bhuiyan L B 1983 J. Chem. Soc. Faraday Trans. II 79 707
DOI:10.1039/F29837900707 [Cited within: 2]

Outhwaite C W Molero M Bhuiyan L B 1991 J. Chem. Soc. Faraday Trans. 87 3227
DOI:10.1039/FT9918703227 [Cited within: 1]

Bari Bhuiyan L Outhwaite C W 2004 Phys. Chem. Chem. Phys. 6 3467
DOI:10.1039/B316098J [Cited within: 2]

Kirkwood J G 1934 J. Chem. Phys. 2 767
DOI:10.1063/1.1749393 [Cited within: 1]

Kirkwood J G Poirier J C 1954 J. Phys. Chem. 58 591
DOI:10.1021/j150518a004 [Cited within: 1]

Outhwaite C W Bhuiyan L B 1982 J. Chem. Soc. Faraday Trans. II 78 775
DOI:10.1039/F29827800775 [Cited within: 1]

Outhwaite C W Lamperski S 2001 Condens. Matter Phys. 4 739
DOI:10.5488/CMP.4.4.739 [Cited within: 1]

Lamperski S Bhuiyan L B 2003 J. Electroanal. Chem. 540 79
DOI:10.1016/S0022-0728(02)01278-0 [Cited within: 1]

Lamperski S Outhwaite C W 2002 Langmuir 18 3423
DOI:10.1021/la011852v [Cited within: 1]

Kjellander R Mitchell D J 1992 Chem. Phys. Lett. 200 76
DOI:10.1016/0009-2614(92)87048-T [Cited within: 1]

Hansen J-P McDonald I R 2006 Theory of Simple Liquids3rd ednHansen J-P McDonald I RNew YorkAcademic291
[Cited within: 1]

Varela L M Perez-Rodriguez M Garcia M Sarmiento F Mosquera V 1998 J. Chem. Phys. 109 1930
DOI:10.1063/1.476770 [Cited within: 1]

Lovett R Stillinger F H 1968 J. Chem. Phys. 48 3869
DOI:10.1063/1.1669710 [Cited within: 1]

Stillinger F H Lovett R 1968 J. Chem. Phys. 49 1991
DOI:10.1063/1.1670358

Stillinger F H Lovett R 1968 J. Chem. Phys. 48 3858
DOI:10.1063/1.1669709 [Cited within: 1]

闂傚倸鍊搁崐鐑芥嚄閸撲礁鍨濇い鏍仜缁€澶嬩繆閵堝懏鍣圭紒鐘靛█閺岀喖骞戦幇闈涙闂佸憡淇洪~澶愬Φ閸曨垰妫橀柛顭戝枓閹稿啴姊洪崨濠庢畷鐎光偓閹间礁绠栨俊銈呮噺閺呮煡骞栫€涙ḿ绠橀柣鈺佹捣缁辨挻鎷呮搴ょ獥闂侀潻缍囩紞浣割嚕婵犳碍鍋勯柣鎾虫捣椤ρ囨⒑閸忚偐銈撮柡鍛箞閹繝宕掗悙绮规嫼缂備礁顑堝▔鏇㈡倿閸ф鐓欓柛鎴欏€栫€氾拷2婵犵數濮烽弫鎼佸磻閻愬搫鍨傞柛顐f礀缁犳壆绱掔€n偓绱╂繛宸簻鎯熼梺鍐叉惈椤戝洨绮欒箛娑欌拺闁革富鍘奸崝瀣亜閵娿儲顥㈢€规洜鏁婚崺鈧い鎺戝閳锋垿鏌涘☉姗堝伐濠殿噯绠戦湁婵犲﹤鎳庢禒杈┾偓瑙勬礃濡炰粙寮幘缁樺亹鐎规洖娲ら獮妤呮⒒娓氣偓濞佳呮崲閸儱纾归柡宓偓濡插牏鎲搁弮鍫濊摕闁挎繂顦悞娲煕閹板吀绨奸柛锝勫嵆濮婅櫣鎷犻垾铏闂佹悶鍎滈崶褎鏆梻鍌欑劍鐎笛呮崲閸屾娲閵堝懐锛涢梺鍦劋椤ㄥ棝鍩涢幋锔界厱婵犻潧妫楅鈺呮煃瑜滈崜娆戠礊婵犲洤绠栭梺鍨儐缂嶅洭鏌嶉崫鍕簽婵炶偐鍠栧铏规崉閵娿儲鐝㈤梺鐟板殩閹凤拷
婵犵數濮烽弫鍛婃叏娴兼潙鍨傜憸鐗堝笚閸嬪鏌曡箛瀣偓鏇㈡倷婵犲嫭鍠愮€广儱妫欓崣蹇涙煏閸繍妲归柍閿嬪灴閺屾稑鈽夊鍫濅紣缂備焦顨嗙敮妤佺┍婵犲浂鏁冮柨婵嗘处閸掓稑顪冮妶鍐ㄧ仾婵☆偄鍟幈銊╁焵椤掑嫭鐓忛柛顐g箖閿涘秵淇婇銏狀伃闁哄矉绲鹃幆鏃堫敍濠婂憛锝夋⒑閸濄儱校闁绘濮撮悾鐑藉閵堝懐顔掑銈嗘⒒閺咁偊宕㈤幖浣光拺闁告稑锕ョ粈瀣箾娴e啿娲﹂崐鍫曟煥濠靛棙顥撳ù婊勭矒閺岀喓鈧稒岣跨粻鏍ь熆鐠哄搫顏紒杈ㄥ笧閳ь剨缍嗘禍璺何熼埀顒勬⒑缁洘鏉归柛瀣尭椤啴濡堕崱妤€娼戦梺绋款儐閹瑰洭寮诲鍥ㄥ珰闁哄被鍎卞鏉库攽閿熺姷鐣哄ù婊冪埣瀵顓奸崼顐n€囬梻浣告啞閹稿鎮烽埡浣烘殾妞ゆ牗绋戦閬嶆倵濞戞顏呯椤栨埃鏀介柣鎰级閳绘洖霉濠婂嫮绠炵€殿喗鐓¢、妤呭礋椤掆偓閳ь剙鐖奸弻锝夊箛椤旇姤姣勯梺纭呮閸婂潡寮诲☉銏犖ч柛銉仢閵忋倖顥嗗璺侯儑缁♀偓婵犵數濮撮崐鎼佸汲閿濆棎浜滈幖娣焺濞堟洟鏌曢崶褍顏柛鈺冨仱椤㈡﹢鎮欏顔荤棯濠电姵顔栭崹閬嶅箰閹惰棄钃熼柨鐔哄Т閻愬﹪鏌嶆潪鎵妽闁诲繋绶氬娲川婵犲嫭鍠涢梺绋款儐閹瑰洤顫忕紒妯诲闁告縿鍎虫婵犵數鍋橀崠鐘诲幢閹邦亝鐫忛梻浣虹帛閸旀寮崫銉т笉闁哄啫鐗婇悡娆撴煙椤栧棗鑻▓鍫曟⒑瀹曞洨甯涙慨濠傜秺楠炲牓濡搁妷搴e枔閹风娀骞撻幒婵囨祰闂傚倷鐒﹂幃鍫曞磹瑜忕划濠氬箻鐠囪尪鎽曢梺缁樻濞咃綁鎯屽▎鎾寸厵缂佸鐏濋銏ゆ煙椤旂晫鎳囨慨濠勫劋鐎电厧鈻庨幋鐘樻粎绱撴担鍝勑i柣妤佹礋椤㈡岸鏁愭径妯绘櫇闂佸啿鐏堥弲婊堟倵婵犳碍鈷戠憸鐗堝笒娴滀即鏌涘Ο鍝勨挃缂侇喗鐟╁畷鐔碱敍濞戞帗瀚奸梻浣告贡鏋繛瀵稿厴閸┿儲寰勯幇顓犲幐闂佸壊鍋掗崑鍕櫠鐎电硶鍋撶憴鍕缂傚秴锕ユ穱濠傤潰瀹€濠冃┑鐘愁問閸ㄤ即濡堕幖浣歌摕婵炴垶菤濡插牊鎱ㄥΔ鈧悧濠囧极閸撗呯=濞达絽鎼牎闁汇埄鍨抽崑銈夊春閳ь剚銇勯幒鍡椾壕闂佽绻戦懝楣冣€﹂崹顕呮建闁逞屽墴楠炲啳顦圭€规洖宕湁闁哄瀵ч崰妯尖偓瑙勬礈鏋摶鏍归敐澶嬫珳闁汇儺浜缁樻媴娓氼垱鏁梺瑙勬た娴滎亜顫忔禒瀣妞ゆ牗绋掑▍鏍⒑閸濆嫮鈻夐柛妯圭矙閹ょ疀濞戞瑧鍘遍梺鏂ユ櫅閸燁垳绮堥埀顒€顪冮妶蹇曞矝闁哄棙绔糴婵犵數濮烽弫鍛婃叏娴兼潙鍨傞柛锔诲幘缁€濠傗攽閻樺弶鎼愰柣鎺戠仛閵囧嫰骞掑鍫濆帯闂佹剚鍨卞ú鐔煎蓟閺囥垹骞㈡俊銈傚亾闁哄棴缍侀弻锛勪沪閸撗勫垱濡ょ姷鍋炵敮锟犵嵁鐎n喗鍊婚柛鈩冿供濡冣攽閿涘嫬浜奸柛濠冪墱閺侇噣鎮欓崫鍕崶闂佸綊鍋婇崰姘舵儗濞嗗繆鏀介柣妯哄级婢跺嫰鏌涚€n偄濮嶉柡宀嬬秮婵偓闁靛繆鍓濆В鍕煛娴e摜澧︽慨濠勭帛閹峰懐绮欓幐搴♀偓顖氣攽閻橆喖鐏柨鏇樺灩閻g兘顢涘☉姗嗗殼闁诲孩绋掗敋濞存粠鍨跺娲川婵犲嫮鐣垫繝娈垮灥妞存悂骞嗛弮鍫濐潊闁挎稑瀚倴濠碉紕鍋戦崐鏍礉濡ゅ懎绐楅幖娣灮椤╂彃螖閿濆懎鏆為柣鎾寸洴閺屾盯濡烽敐鍛瀴闂佹眹鍔嶉崹鍧楀蓟閿濆鍋勯柛娆忣槹閻濇棃姊虹€圭姵顥夋い锔炬暬閻涱喖螣閼测晝顦╅梺缁樏畷顒勵敆閵忊€茬箚闁绘劦浜滈埀顒佺墪鐓ゆ繝闈涙閺嬪秹鏌¢崶鈺佷憾缂傚倹宀搁悡顐﹀炊閵娧€妲堥悗鐟版啞缁诲啴濡甸崟顖氱婵°倐鍋撻柛鐕佸灦椤㈡瑩鏁撻敓锟�20濠电姴鐥夐弶搴撳亾濡や焦鍙忛柣鎴f绾惧鏌eΟ娆惧殭缂佺姴鐏氶妵鍕疀閹炬惌妫″銈庡亝濞叉ḿ鎹㈠┑瀣棃婵炴垵宕崜鎵磽娴e搫校闁搞劌娼″濠氬Χ閸℃ê寮块梺褰掑亰閸忔﹢宕戦幘婢勬棃鍩€椤掑嫬鐓濋柡鍐ㄧ墕椤懘鏌eΟ鐑橆棤闁硅櫕鎹囬妶顏呭閺夋垹顦ㄩ梺鍐叉惈閿曘儵鏁嶉崨顖滅=闁稿本鐟чˇ锔姐亜閿旇鐏︽い銏″哺椤㈡﹢濮€閻橀潧濮︽俊鐐€栧濠氬磻閹惧绡€闁逞屽墴閺屽棗顓奸崨顖ょ幢闂備胶绮濠氬储瑜斿鍛婄瑹閳ь剟寮婚弴銏犻唶婵犲灚鍔栨晥闂備胶枪妤犲摜绮旇ぐ鎺戣摕婵炴垯鍨归崡鎶芥煏婵炲灝鍔氭い顐熸櫊濮婄儤瀵煎▎鎴犳殸缂傚倸绉撮敃顏堢嵁閸愩剮鏃堝礃閳轰焦鐎梻浣告啞濞诧箓宕f惔銊ユ辈闁跨喓濮甸埛鎴︽煕濠靛棗顏い銉﹀灴閺屾稓鈧綆鍋呭畷灞炬叏婵犲啯銇濈€规洦鍋婂畷鐔煎垂椤愬诞鍥ㄢ拺闁告稑锕ラ埛鎰版煟濡ゅ啫鈻堟鐐插暣閺佹捇鎮╅搹顐g彨闂備礁鎲″ú锕傚礈濞嗘挻鍋熷ù鐓庣摠閳锋垿姊婚崼鐔恒€掔紒鐘冲哺閺屾盯骞樼€靛摜鐤勯梺璇″枓閳ь剚鏋奸弸搴ㄦ煙闁箑鏋ゆい鏃€娲樼换婵嬪閿濆棛銆愬銈嗗灥濡稓鍒掗崼銉ョ劦妞ゆ帒瀚崐鍨箾閸繄浠㈡繛鍛Ч閺岋繝鍩€椤掑嫬纭€闁绘垵妫楀▓顐︽⒑閸涘﹥澶勯柛瀣浮瀹曘儳鈧綆鍠楅悡鏇㈡煛閸ャ儱濡兼鐐瓷戞穱濠囧矗婢跺﹦浼屽┑顔硷攻濡炶棄鐣烽锕€绀嬫い鎰枎娴滄儳霉閻樺樊鍎滅紓宥嗙墪椤法鎹勯悜妯绘嫳闂佺ǹ绻戠划鎾诲蓟濞戙埄鏁冮柨婵嗘椤︺劑姊洪崫鍕闁告挾鍠栭獮鍐潨閳ь剟骞冨▎鎾搭棃婵炴垶顨呴ˉ姘辩磽閸屾瑨鍏屽┑顔炬暩閺侇噣鍨鹃幇浣圭稁婵犵數濮甸懝楣冩倷婵犲洦鐓ユ繝闈涙閸gǹ顭跨憴鍕婵﹥妞介幊锟犲Χ閸涱喚鈧儳鈹戦悙鎻掔骇闁搞劌娼¢獮濠偽旈崘鈺佺/闁荤偞绋堥崜婵嬫倶娓氣偓濮婅櫣娑甸崨顔兼锭闂傚倸瀚€氭澘鐣烽弴銏犵闁挎棁妫勯埀顒傛暬閺屻劌鈹戦崱娑扁偓妤侇殽閻愮榿缂氱紒杈ㄥ浮閹晛鐣烽崶褉鎷伴梻浣告惈婢跺洭宕滃┑鍡╁殫闁告洦鍋€濡插牊绻涢崱妤佺濞寸》鎷�
相关话题/brief review continuous