# Modified Poisson-Nernst-Planck model with accurate Coulomb correlation in variable media

###### Abstract

We derive a set of modified Poisson-Nernst-Planck (mPNP) equations for ion transport from the variation of the free energy functional which includes the many-body Coulomb correlation in media of variable dielectric coefficient. The correlation effects are considered through the Debye charging process in which the self energy of an ion is governed by the generalized Debye-Hückel equation. We develop the asymptotic expansions of the self energy taking the ion radius as the small parameter such that the multiscale model can be solved efficiently by numerical methods. It is shown that the variations of the energy functional give the self-energy-modified PNP equations which satisfy a proper weak formulation. We present numerical results from different asymptotic expansions with a semi-implicit conservative numerical method and investigate the effect of the Coulomb correlation.

harge transport; ion correlation; continuum theory for electrostatics; asymptotic expansion; Green’s function; Poisson-Nernst-Planck equations

82C21; 82D15; 35Q92; 49S05

## 1 Introduction

The understanding of the transport of charged particles in micro/nanoscale systems is of importance in many natural phenomena and has wide engineering applications [1]. The Poisson-Nernst-Planck (PNP) theory, which is composed of the Nernst-Planck equations for the ionic flux under the concentration gradient and electrical field, and the Poisson equation for the electrical potential, is a classical continuum model to describe the ionic transport in solution and has been shown to be successful in lots of applications including nanofluidics and microfluidics [2, 3], semi-conductor devices [4] and transmembrane ion channels [5].

The traditional PNP model cannot make correct prediction for systems with divalent ions or with a dielectric boundary because it treats ions as point charges without short-range size effect and ignores the long-range Coulomb correlation. The improvement of the mean-field approximation has attracted wide attention in recent decades [6, 7, 8, 9, 10, 11, 12, 13, 14]. The steric effect between particles can be considered by introducing the excess free energy functional given by solvent entropy [10], Lennard-Jones kernel [9], Carnahan-Starling local density approximation [15, 16] or modified fundamental measure theory [17, 18], see reviews [19, 13] for discussions of this issue. For the long-ranged Coulomb correlation, Bazant et al. [20] proposed a fourth-order Poisson-Boltzmann (PB) model by introducing a parameter of correlation length, which has been further developed recently [21, 22, 23]. Another promising routine is by introducing a self-energy correction to the potential of mean force using the generalized Debye-Hückel equation, which can be derived under the point-charge approximation (PCA) by the Debye closure of the BBGKY hierarchy [24], random phase approximation in density functional theory [25] or the Gaussian variational field theory [26, 27, 28] for equilibrium systems. The PDEs under the PCA are unstable [29] for strong-coupling systems because cations and anions may collapse when the interaction is strong. To go beyond PCA, Wang [30] used the variational field theory to derive a generalized Debye-Hückel equation for finite-sized ions with Gaussian charge distribution, which first revealed the importance of the ionic self energy for mobile ions with finite size and obtained the Born solvation energy in variable dielectric media by performing expansion using ionic radius as small parameter, and provides the theoretical framework for accounting for these many-body effects. It should be noted that the dielectric permittivity in electrolytes could depend on the ionic concentration as well as the electrical field [31, 32, 19, 33, 34, 35], thus the Born solvation energy should play important role due to strong variations in the ion concentration or the electric field. Accounting for the Born energy and ionic correlation in a continuum model is essential to explain many phenomenon such as the nonmonotonic concentration dependence of the mean activity coefficient of simple electrolytes [36] and quantitative investigation of single ion activities [37]. The extension of these models from equilibrium to time-dependent dynamics with the Fick’s law is straightforward, and results in the modified PNP equations.

In this work, we derive the system free energy through the Debye charging process [38], which accounts for the Coulomb correlation of ions with finite size and the dielectric effect through a self energy term. The self energy is described by the generalized Debye-Hückel equation with variable coefficients. When finite ion size is considered, the equation is high-dimensional as the local dielectric coefficient is modified by the located ion, i.e., six dimensions for real 3 dimensional systems. Furthermore, the difference between ionic size and the Debye screening length leads to a multi-scale problem. To overcome these difficulties, asymptotic expansions of the self energy by taking the ionic size as a small parameter. The resulting modified PB/PNP equations or similar versions had been used by different authors to investigate important physical features such as the double layer structure and like-charge attraction [39, 40, 41, 42]. With the asymptotic self-energy expression, we find the connection between the chemical potential and self energy, so that the modified PNP equation can preserve a simple form and satisfy a proper energy dissipation law thus a weak formuation. We compare the results by different kind of approximations, and demonstrate the new formulations can predict correct physics in more accurate.

On the other hand, numerical solutions of the PNP and the modified PNP equations can be very challenging, and have attracted much interest recently [43, 44]. The modified PNP equations subject to suitable boundary conditions make new difficulties in numerical solutions, since the correlation function may lead to the local condensation of ions and the non-monotonic ionic distributions. In Section 3, we develop a semi-implicit numerical method to solve the modified PNP equations with the mass conservation. We perform numerical example to show that the method can capture correct physical features.

## 2 Coulomb correlation in variable media

To take into account the correlation contribution in the PNP framework, we start from the free energy formulation. The modified PNP equations are derived through the energetic variational approach [9, 45] such that the energy dissipation law and a proper weak formulation will be automatically fulfilled.

Let be the electrolyte region and be its boundary. Suppose there are ionic species in the electrolyte with being the valence of species and each ion is modeled as a hard sphere with radius and a point charge in the center. In the mean-field theory, the free energy functional is given by,

(1) |

where is the concentration distribution of the th ion species, is the Boltzmann constant, is the temperature, is the space dependent dielectric permittivity, and the mean electrical potential satisfies the Poisson’s equation,

(2) |

with proper boundary conditions. Since the dielectric environment is inhomogeneous, we should consider the contribution from the Born energy of each ion,

(3) |

with defined below in Eq. (6). Therefore the finiteness of the ionic radius is essential for variable media as the zero-radius limit leads to the divergence of the Born solvation energy. Here is not assumed to be concentration dependent nor field dependent for simplicity, or we need extra terms to describe the polarization.

### 2.1 Debye charging process

The Debye charging process considers a hypothetical process, in which all ions are charged with the same rate from zero charge to its full charge. The work done during the process is equal to the free energy difference and thus can be used to calculate ionic chemical potential or activity coefficients [38]. When a source charge is located at , the dielectric permittivity is locally dependent on the source location due to the finite size of the ion,

(4) |

We choose the reference system to be the ideal gas, i.e., its particle-particle interaction is , and increase the valences of all particles up to full charges: At stage with , the valence of the th species becomes , and the pairwise potential is thus,

(5) |

where is governed by,

(6) |

When , it is the ideal gas reference system, and when , it becomes the real electrolyte system of interest.

We introduce the standard perturbation theory for many-body systems [46]. Let be the correlation free energy as a functional of single particle density. Then it can be expressed as,

(7) |

where is the total correlation function, with being the two-particle density distribution function, which can be viewed as the joint distribution of and . In the mean-field theory and , meaning no correlation is included. In our model, we are motivated to go beyond mean-field theory and take into account the electrostatic correlation between ions. We describe the total free energy as,

(8) | |||||

The second term is the correction to the classical mean field theory, where the self-energy includes the Born solvation energy (3) and electrostatic correlation (7),

(9) |

### 2.2 Self energy

The self energy (9) is treated at the Debye-Hückel level, which means only the test ion is considered with finite radius while the surrounding mobile ions are still treated as point charges. Let be the spherical domain of a test ion located at , and is the indicator function of its complementary set , i.e., for and 1 elsewhere. Assume adding an hard sphere at location will not change the surrounding ionic distribution except the interior domain of the ion sphere is inaccessible, i.e. the short-range hard core correlation is ignored. Then replace the hard sphere with the th ion, as a consequence, the electrostatic potential will change by , which satisfies,

(10) |

where the first term on the right hand side is the charge fluctuation due to the insert of the ion and the second term represents the point charge within the ion. We use the Loeb closure [47] in domain ,

(11) |

where is the inverse thermal energy. Then we obtain the generalized Debye-Hückel (GDH) equation [48, 40] with ,

(12) |

where is the ionic strength defined by,

(13) |

Notice that the pairwise potential between the th fixed test ion and the th free ion is . The self energy of the th ion defined in Eq. (9) is expressed as,

(14) |

Now we have the total free energy as a functional of the ionic concentrations, and the chemical potential can be computed through the functional derivatives. In literature, there is another process called the Güntelberg charging process [49], where only one central ion goes through the hypothetical charging, while other ions remain fully charged. Similar derivation can also be applied to the Güntelberg charging process, resulting in the excess chemical potential,

(15) |

However, these two charging processes do not necessary give the same result of chemical potential [50] and unfortunately, the chemical potential given by the free energy functional (8) does not have the simple form as (15).

### 2.3 Asymptotic expansions

In order to develop a PDE model that can be solved efficiently, we apply the asymptotic method for (12) to obtain a simplified equation by considering the ionic radius to be a small parameter in comparison with the inverse Debye length. We first introduce the PCA which together with the Born energy is the zeroth order asymptotics of the GDH equation, and then discuss higher order extensions.

#### 2.3.1 Point charge approximation

At the limit of , ions are point charges and the GDH equation becomes,

(16) |

Since the Born energy becomes singular at the zero radius, we shall consider it later on and define the self energy by,

(17) |

The Green’s function can be viewed as the inverse of the modified Helmholtz operator, which is still hard to solve by numerical methods. Since we are only interested in the diagonal components, i.e., the self Green’s function , the selected inversion algorithm [51] can be employed for the purpose [29]. Although the self energy is well-defined, however it ignores the excluded-volume effect and would lead to the instability for strong correlated systems due to the ionic collapse. Furthermore, since it does not include the Born energy of an ion, an appropriate treatment of it should introduce the ionic radius again.

#### 2.3.2 Higher-order approximations

The Green’s identities are employed to expand the asymptotic expansion into higher orders. Let us rewrite the GDH equation (12) into an integral form,

(18) | |||||

where and are divergence and gradient operators with respect to . Similarly,

(19) | |||||

We shall use the fact that outside and inside , and notice that is the indicator function of . Subtracting (19) from (18) and use Green’s identities, we have,

(21) | |||||

where is the outer normal derivatives. We use the interface conditions on , namely, functions , , and are all continuous across the interface, and we then obtain,

(22) | |||||

where is the inverse of the local Debye-length function, . Define the nonsingular parts of and through,

(23) |

Thus, in whole space of , satisfies,

(24) |

Assume and are smooth functions, then and are smooth too. Since we are only interested in the diagonal entries of , by substituting (23) into (22) and taking the limit , we obtain,

(25) |

where , and are the integrals given in (22) after taking the limit , i.e.,

(26) |

Now we can obtain a higher-order asymptotic expansion by calculating the integrals within the ball and the integral on its boundary. For the leading order, we can use the singular parts of and . This yields,

(27) |

Clearly, , , and has the leading asymptotics of order ,

(28) |

which dominates the right hand side of Eq. (25). Since , we can easily obtain a zeroth-order approximation of the self energy,

(29) |

This is the Born energy correction to the PCA.

In order to compute the coefficients of terms, we need to evaluate the integrals up to one more order. We have,

(30) |

where

(31) |

Substituting these terms into Eq. (25), we arrive at the first-order asymptotics,

(32) |

In the expression, the first term is the PCA correlation energy, the second term is the Born energy due to the dielectric variation, and the two terms in are corrections to the PCA due to the finite ion size and to the Born energy, respectively. When the dielectric coefficient is highly inhomogeneous, the corrections become important. In particular, the correction to the Born energy is often ignored in the continuum theory and Monte Carlo simulations of charged systems in variable media [52, 35], and its effect on the ionic structure is less investigated.

The higher-order asymptotics can be obtained by further expanding the three integrals. This will lead to a tedious calculation. We only compute the correction term from the ionic correlation and ignore the contribution from inhomogeneous dielectric media in the second-order term, which gives an additional correction,

(33) |

to Eq. (32). We then obtain a second-order asymptotics, , for systems with slow dielectric variations. It can be seen that the smallness of is in comparison with the local Debye length, i.e., .

#### 2.3.3 Accuracy validation of the asymptotics

To validate the accuracy of these asymptotic approximations, we study two special cases and present numerical calculations for a toy example with more complicated setup.

##### Case 1

In the original Debye-Hückel theory, the self energy is calculated in the bulk region where and the ionic strength are both constants. In this case, the solution of (12) can be computed analytically,

(34) |

In this case, the PCA Green’s function is the screened Coulomb potential and thus we can easily find the PCA self energy The exact self energy and its asymptotic expansion are,

(35) | |||||

This analytical result is consistent with our asymptotic approximations. The comparison results are shown in Fig. 1(a).

##### Case 2

In this case, we fix and consider the case that and are only functions of radial distance and satisfy the relation,

(36) |

where is a constant. We can solve the solution of PCA Green’s function with the source at the origin ,

(37) |

The PCA self energy is The solution of the GDH equation (12) is,

(38) |

which yields the self energy,

(39) |

Since is smooth and at least the second-order derivative exists, we can expand it to obtain, We then have,

(40) |

which is consistent with the first order approximation (32).

##### Example 1

In this example, the dielectric permittivity is piecewise constant: inside the ion, , when , , and otherwise . we choose to mimic the ionic distribution near a charged ion immersed in an electrolyte. We calculate the self energy, with different asymptotic approximations, including the 0th, the 1st and the 2nd corrections. The exact data are prepared by numerical solution of the 3D partial differential equation (12) by the finite-difference method [53] and the reaction potential is computed with a mesh size with grid points. In order to check how accurate this mesh size is, we calculate the numerical solution of Case 1 which is in agreement with the exact solution (solid line) from Fig. 1(a). The results in Fig. 1(b) show that the 2nd correction asymptotics with a varying ionic strength and dielectric coefficient match the exact data well, and the zeroth order asymptotic solution significantly deviates from the exact data. It is also shown that the 1st order correction becomes less accurate with the increase of the ionic size, and when is large the 1st order correction is not enough to provide high accuracy.

### 2.4 Self-energy-modified PNP model

The self energy defined through the point charge approximation is a functional of ionic concentration. In order to obtain the transport equations for the ions, we need to compute the functional derivative of free energy with respect to ionic concentration. First we will present some useful lemmas.

The functional derivative of the PCA self energy (17) with respect to ionic concentration is given by,

(41) |

Suppose that is the perturbed concentration, where is small, the corresponding ionic strength is denoted as , and the PCA self energy after perturbation is described by,

(42) |

Dropping the higher order terms, the difference for the PCA Green’s function equation between (16) and (42) is,

(43) |

The solution can be then expressed by the Green’s identity,

(44) |

Thus the functional derivative leads to Eq. (41).

From this lemma, we can easily see the PCA self energy satisfies the following reciprocal relation,

(45) |

The derivative of the PCA self energy (17) with respect to the coupling parameter is given by,

(46) |

Consider the system with a small perturbation in , . In the new system, the generalized Debye-Hückel equation and the self energy become,

(47) |

Then the difference between Eqs. (47) and (16) leads us to,

(48) |

We then have,

(49) |

Taking the derivative with respect to , we get the result of expression (46).

The Debye charging process of the PCA self energy gives the correlation free energy whose corresponding chemical potential is for the th species.

The correlation energy,

(50) |

The corresponding excess chemical potential is given by the functional derivative,

(51) |

The second term in Eq. (51) is,

(52) | |||||

where the first two equalities use Lemmas 2.4 and 2.4, respectively, and the last equality is simply the result of the integration by parts. The chemical potential can be then written as,

(53) |

which ends the proof of this Lemma.

Lemma 2.4 implies that the chemical potential is exactly the self energy at the full charging , and thus the Debye charging process is consistent with Güntelberg charging process in the PCA situation. With this lemma, we can easily obtain the simplified expression for the correlation energy in the higher order approximations.

Let

(54) |

be the correlation component of the free energy with , 1 or 2 to represent the order of asymptotics. We have:

(1) The zeroth-order and first-order excess chemical potentials are,

(55) |

(2) The second-order asymptotics does not have the same conclusion, but it holds the reciprocal relation,

(56) |

which implies there exists a free energy functional whose variation is the self energy at the full charging . We will redefine .

(1) Since the Born energy is independent of the ionic concentration, we have,

(57) |

Noticing that , we can obtain,

(58) |

Using Lemma 2.4 for the PCA term, we find that the excess chemical potential equals the self energy of an ion up to the first-order asymptotic expansion. We complete the proof of Eq. (55).

(2) For the second-order approximation, following the same procedure as (51), we can easily see that the increment in chemical potential does not equal the increment in the self energy at full charing,

(59) |

The PCA energy satisfies the reciprocal relation. From the proof of Lemma (2.1), we can straightforwardly obtain that,

(60) |

And now for the second-order correction,

(61) | |||||