# Self-energy-modified Poisson-Nernst-Planck equations: WKB approximation

and finite-difference approaches

###### Abstract

We propose a modified Poisson-Nernst-Planck (PNP) model to investigate charge transport in electrolytes of inhomogeneous dielectric environment. The model includes the ionic polarization due to the dielectric inhomogeneity and the ion-ion correlation. This is achieved by the self energy of test ions through solving a generalized Debye-Hückel (DH) equation. We develop numerical methods for the system composed of the PNP and DH equations. Particularly, towards the numerical challenge of solving the high-dimensional DH equation, we developed an analytical WKB approximation and a numerical approach based on the selective inversion of sparse matrices. The model and numerical methods are validated by simulating the charge diffusion in electrolytes between two electrodes, for which effects of dielectrics and correlation are investigated by comparing the results with the prediction by the classical PNP theory. We find that, at the length scale of the interface separation comparable to the Bjerrum length, the results of the modified equations are significantly different from the classical PNP predictions mostly due to the dielectric effect. It is also shown that when the ion self energy is in weak or mediate strength, the WKB approximation presents a high accuracy, compared to precise finite-difference results.

###### pacs:

82.45.Un, 04.25.Nx, 82.60.Lf, 02.70.Bf^{†}

^{†}preprint: Preprint

## I Introduction

The charge transport in fluids under confinements or around objects of nanometer length scale has been of growing interest in a lot of physical and biological systems Schoch et al. (2008); Daiguji (2010); for example, in the study of colloidal separation and self-assembly, nanoparticles at liquid-liquid interfaces, electrochemical energy devices and membrane ionic channels. When the length scale of confinements is comparable with the characteristic lengths of electrolytes (e.g., the Bjerrum length of the solvent , and the Debye length of the electrolyte ), complex electrostatic phenomena such as the charge inversion and the like-charge attraction have been often observed experimentally. These phenomena are beyond mean-field theoretical explanations and motivate many challenges to computational and modelling communities French et al. (2010).

In the vicinity of a charged surface, counterions are attracted, forming a screening region called electric double layer. The structure of the double layer plays a key role in nanoscale interface sciences Levin (2002); Walker et al. (2011); French et al. (2010); Messina (2009); Wan et al. (2014). The ions in the external diffuse layer is Coulombic, where the electric potential decays exponentially with the characteristic Debye length. The ions in the internal Stern layer are condensed, where the electric potential does not obey the exponential rule due to the tangential ion correlation and other short-range interactions. If the medium contacting the electrolyte is a dielectric of low permittivity, the dielectric mismatch leads to an induced-charge potential, repelling the counterions and thus creating a depletion zone near the surface. This effect could play an important role in many phenomena of nanoscale systems, and is also attracting much attention in theoretical and simulation study (see Messina (2002); Hatlo and Lue (2008); Jho et al. (2008); Wang and Ma (2010); Bakhshandeh et al. (2011); Diehl et al. (2012); Gan et al. (2012); Xu (2013); Wang and Wang (2013); Zwanikken and de la Cruz (2013) to mention a few recent literature).

In the mean-field description of ion structure and transport in the double layer, the Poisson-Boltzmann (PB) and Poisson-Nernst-Planck (PNP) theories have been considered. The PNP is beyond the PB because the PNP describes the charge dynamics and thus has been widely used for studying ion transport in nanopores and nanochannels; for example, open ion channels in cell membranes Chen and Eisenberg (1993); Eisenberg (1996, 1999); Hollerbach et al. (2001). The system of equations has also been widely applied in continuum theories of narrow channels by coupling with the Stokes equations to model nanofluidics Daiguji et al. (2005); Daiguji (2010). The PNP equations are also known as “drift-diffusion equations” and have much use in semiconductors Markowich et al. (1990). Asymptotic analysis on the PNP equations has been made to understand properties of charge diffusion of electrochemical energy systems Bazant et al. (2004); Squires and Bazant (2004); Singer and Norbury (2009). From the aspect of numerical analysis, various numerical methods have been proposed Kurnikova et al. (1999); Cardenas et al. (2000); Zheng et al. (2011); Wei et al. (2012); Flavell et al. (2014); Dyrka et al. (2013) to solve the nonlinear equations more accurately and efficiently and to capture specific dynamical behaviors of the non-equilibrium system.

The classical PNP theory neglects the ion-ion correlation, excluded volume and image charge effects in the double layer, significantly changing the local structure and the far-field properties of electrolytes. When large surface charges or applied potentials are involved, the nonmonotonic differential capacitance of electrodes can not be predicted without these effects Lockett et al. (2008). Analysis of current-concentration curves Mamonov et al. (2003) also showed significantly different trends between the experimental data Busath et al. (1998) and that predicted by the PNP theory. Moreover, it has been found that the mean-field theory fails in explaining, even qualitatively, phenomena of the long range attractions between two surfaces in electrolyte Kékicheff and Spalla (1995). The comparison between the PNP theory and particle-based Brownian dynamics simulations Corry et al. (2000) for cylindrical pores of varying sizes demonstrated the invalidity of the mean-field theories when the cylindrical radius is less than , evidencing that the aforementioned effects should be accounted for in the continuum modeling.

Different versions of modified PB/PNP theory have been proposed by including such as steric effects Storey et al. (2008); Kilic et al. (2007); Horng et al. (2012); Lu and Zhou (2011); Eisenberg et al. (2010); Liu et al. (2012) or the dielectric self energy Nadler et al. (2003); Corry et al. (2003); Graf et al. (2004); Hyon et al. (2011). A useful approach of remedying the mean-field theory is to replace the mean potential by a better approximation of the potential of mean force, i.e., to include the ignored effects by correcting the mean potential by the self energy of mobile ions. The remedy to modify the mean potential may improve the mean-field theory a lot by including ionic correlations and image charges to capture many-body physics; see, e.g., Luo et al. (2006); Wang and Wang (2013). The shortcoming is that the self energy is very difficult to obtain, either by molecular dynamics simulations or by solving high-dimensional equations from field-theory calculations. In this work, we follow the idea of Gaussian variational field theory Podgornik (1989); Netz and Orland (2000, 2003) which models the self energy as the self-Green’s function. A self-energy modified PNP model is derived (Sec. II) and the Green’s function is described by the generalized Debye-Hückel (DH) equation nonlinearly depending on the mean ion distributions. The resulted system of partial differential equations is then composed of three equations: the Nernst-Planck equations for the dynamics of charges, the Poisson equation for the mean potential, and the generalized DH equation for the Green’s function. The self energy is solved efficiently by some analytical and numerical approximations, and then the modified PNP equations are solved by finite difference algorithm.

In Section III, we discuss the dimensionless formulation of the modified PNP equations, by introducing two ratios of length scales, and , where is the characteristic length of confinements. Then we developed the analytical WKB (Wentzel-Kramers-Brillouin) approximation and the numerical finite-difference approximation for the DH equation. The WKB approach attempts to express the self energy as an explicit formula by asymptotic approximation assuming the parameter is small. The finite-difference numerical method is more expensive, but can be tackled by using selective inversion for sparse symmetric and positive definite matrix. In Section 4, our numerical results show the WKB approximation is in a good agreement with the difference method when is bigger than a few Bjerrum lengths. This is a positive evidence that the use of the WKB approximation gives a satisfactory accuracy in the self-energy calculation, while avoiding expensive high-dimensional numerical calculations.

Our numerical examples place much focus on the dielectric effect of electrodes. This effect is believed to be important in many systems as aforementioned, and may play important role in many-body interaction between colloids Wang and Wang (2013); Zwanikken and de la Cruz (2013). Certainly, electrostatic interaction at solid-liquid interfaces is much more complex than the simplified picture. One important effect is the variable dielectric permittivity of solvent, which has been found to capture many key phenomena observed experimentally Bonthuis et al. (2011); Bonthuis and Netz (2013). The modified PNP model does include the treatment of this effect with varying dielectric constant and excluded volumes of test ions by the DH equation; Eq. (4) below. This treatment will highly increase the computational cost compared with the present model, thus, we remain our model here for a dielectrically homogeneous electrolyte bounded by sharp interfaces.

## Ii Governing equations

The transport of charged particles in an electrolyte is often described by the Nernst-Planck equation (also called the Smoluchowski equation) which states that the time derivative of the ion density function is composed of a diffusion contribution and an advection contribution due to the potential energy of the ions,

(1) |

where is the concentration of the ions of species , is its diffusion constant, where is the Boltzmann constant and is the temperature.

In the mean-field description of electrostatics, the potential energy takes the mean potential energy , and the electric potential is governed by the Poisson equation,

(2) |

where is the valence, is the fixed charge, is the vacuum dielectric permittivity and is the relative dielectric permittivity of the medium. The coupling system between the Nernst-Planck equation and the Poisson equation is called the Poisson-Nernst-Planck (PNP) equations.

The mean-field nature of the PNP ignores the induced-charge effect (or image-charge effect) of dielectric discontinuity and many-body ion correlation, which plays an important role in a lot of electrostatic phenomena. Even at the weak-coupling limit, the induced-charge effect is significant and should be taken into account. In order to include these effects, the potential energy of transported particles can be expressed as the mean potential energy plus a correction,

(3) |

where is the self energy of a unit test ion of the th transported species. Given the mean potential , an accurate description of is to include all many-body interaction with the test ion, i.e., the potential of mean force, which for example can be done by molecular dynamics simulations Luo et al. (2006), or by a reaction-field formulation Corry et al. (2003). In the self-consistent Gaussian field approximation Netz and Orland (2000, 2003); Wang (2010), this quantity can be defined through a Green’s function , described by the generalized Debye-Hückel (DH) equation,

(4) |

and is then the self Green’s function limit,

(5) |

where is the Green’s function in free space to remove the invariable singularity, and is the local ionic strength, which describes the screening effect by the surrounding ions of the test ion, and where the ionic concentrations are determined by the Nernst-Planck equation. The prime in the permittivity function describes the function has been locally modified due to the excluded volume of the test ion,

(6) |

where is the ionic radius, and is the effective dielectric permittivity inside the ionic cavity. We see the ionic size effect is taken into account by assuming that the ionic cavity is inaccessible to mobile ions. Consequently, in Eq. (4), is expressed as,

(7) |

This modification is particularly useful to remove the self-energy singularity when the permittivity of solvent is space-dependent Cherepanov et al. (2003); Bonthuis and Netz (2013), but greatly increases numerical difficulty. To avoid this difficulty, we will not discuss this size effect in this work, and simply take the limit , and leave the algorithm development for a future study.

In equilibrium, the zero ion flux of each species leads to the following equality,

Solving this equation gives an explicit formula for the equilibrium ion density,

where is constant determined by the chemical potential. Substituting the expression into the right side of the Poisson equation, we obtain a modified Poisson-Boltzmann equation,

(8) |

Together with the generalized DH equation (4), it has been studied by Avdeev and Martynov Avdeev and Martynov (1986) using the Debye closure of the BBGKY chain, Netz and Orland Netz and Orland (2000, 2003) using Gaussian variational field theory, and in much recent work (see Buyukdagli et al. (2012); Yaroshchuk (2000); Hatlo and Lue (2008); Lau et al. (2002) to mention a few).

## Iii Charge dynamics in the presence of planar surfaces

In order to understand the self-energy effects for the charge transport, we consider an electrolyte with 1:1 salt between two parallel planar electrodes at (see Fig. 1), which is a simple model of electrochemical systems. We study the case of sharp dielectric permittivity which takes the water permittivity for with a small , and the alternative value for the electrodes outside. The use of a small separation between the dielectric interface and the electrode avoids the self-energy divergence near the boundary. The sharp interfaces remain, providing the induced-charge effect to the mobile ions in the electrolyte.

### iii.1 Dimensionless equations, boundary and initial conditions

Let be the Bjerrum length at the water solvent, at which distance the interaction energy of two unit point charges is . We take as the reference length scale, as a typical diffusion constant, and as a typical ion condensation. We assume ions have uniform diffusion constant and define the another length scale . We note, when is the bulk ion concentration for symmetric monovalent salt, is the Debye screening length.

The following dimensionless parameters and variables are defined: , and , , , , , . We drop the tildes of all new variables. Then the modified PNP and generalized DH equations are,

(9) | |||

(10) | |||

(11) | |||

(12) | |||

(13) |

where and are two length-scale parameters. We see and describe the effective diffusion constant of the Nernst-Planck equation, and dielectric permittivity of the Poisson and generalized DH equations, and represents the charge of the test ion, i.e., the strength of the self energy.

The Nernst-Planck and the Poisson equations are defined on a finite interval , while the generalized DH equation is defined on the whole three-dimensional domain with implied interface conditions at . We assume completely blocking electrodes Bazant et al. (2004) imposed with a time-varying external potential, in dimensionless unit, on the electrodes. Then, at , the ionic fluxes should vanish and the electric potential is a fixed time function. This boundary condition leads to a global zero flux at the steady state, i.e., the equilibrium solution. The current-voltage relations for channel problems correspond to different boundary conditions and will be explored in a future publication under the self-energy modified PNP model. Concerning the generalized DH equation, the Green’s function should be solved in an infinite domain, i.e., the decaying boundary condition. Therefore, we have the following boundary conditions for concentrations, potential and the Green’s function,

(14) | |||

(15) | |||

(16) |

In this work, we take to use the Dirichlet condition for the potential. The Robin boundary condition, i.e., , is often used to account for the dielectric electrodes in literature Bazant et al. (2004); Flavell et al. (2014) by considering the tight counterion adsorption at the Stern layer. In the Green’s function equation, the dielectric permittivity is discontinuous at interfaces at . We find as a result of using a positive , the higher-order singularity of the self energy can be removed when . The dimensionless dielectric permittivity is then for , and elsewhere.

Initially, it is assumed that the ions are in equilibrium without applying the external potential which starts to function at . Therefore, we can set constant initial ionic concentrations, , under electric-neutral constraint . For monovalent salt we are studying, for and 2. In addition, we assume there is no fixed charge in the system, though fixed charges play important roles in many biological and nanofluidic devices.

Now we have completed initial and boundary conditions for the PNP model. The system properties are determined by the ionic species, the applied external potential, and the ratios of length scales and .

### iii.2 Discretization of the PNP equations

The Nernst-Planck and the Poisson equations are approximated by finite-difference discretizations. Let and be the time and space steps, and denote the th time and th space grids by and , respectively. The second-order semi-implicit time-stepping scheme is employed for the Nernst-Planck equation,

(17) |

where and Here the scheme of diffusion term is a central difference, and that of the advection term is an extrapolation for approximating . By this discretization, we gain the second-order accuracy in time and the benefit of avoiding nonlinear iterations thanks to the use of the explicit approximation for the nonlinear advection coefficient. The stability condition mainly depends on the implicitly-discretized diffusion term, and thus the grid sizes can be . Since it is a three-point scheme in time, we should use the backward scheme for at the initial step .

Concerning the spatial discretization, because the dielectric function is constant in interval and the ionic concentrations are smooth, we could rewrite the advection term as . Standard three-point central differences are used for the first and the second space derivatives of and for the internal points. In discretizing the boundary conditions, two ghost points outside the boundaries are introduced to obtain second-order central approximations. For the Poisson equation with given , the mean potential at th time step can be simply obtained by central discretization with the boundary condition.

### iii.3 Solution of the generalized DH equation

Different from the Nernst-Planck and the Poisson equations, the discretization for the generalized DH equation is not trivial. We rewrite the generalized DH equation as a space-dependent-coefficient form, which is required to solve in each time step,

(18) |

where the generalized inverse Debye length, between two interfaces and zero otherwise. Here we set without loss of generality. Solving this equation gives the Green’s function at time . This equation is more difficult due to higher dimensions of the Green’s function.

We will introduce two approaches. The WKB approximation is often adopted to find approximate Green’s function in the presence of electrolytes and interfaces due to its analytical nature Buyukdagli et al. (2012); Wang and Wang (2013), which avoids numerical solution of the high-dimensional problem. We propose a simple WKB expression by improving the idea of Buff and Stillinger Buff and Stillinger (1963) which has been used recently for studying the double layer interaction by Wang and Wang Wang and Wang (2013). This approach results in an explicit formula of the self energy with clear physical significance. In the second approach we propose a numerical approximation with finite difference, which could reach any desired accuracy by varying the grid sizes, and so we limit the use of WKB to a pure analytical formulation.

#### iii.3.1 WKB approximation

In the WKB approximation introduced by Buff and Stillinger Buff and Stillinger (1963) for one-interface problem, the Green’s function is first exactly found in the case of being zero, then the approximate Green’s function takes the screened Coulomb potential for each image charges using the local ion concentration for the inverse screening length. In the presence of two interfaces, the Green’s function of the salt-free solution is a series of image charges by the reflection between two interfaces, as illustrated in Fig. 1. For finite , by the WKB, we have the Green’s function as,

(19) |

where is the distance between and the th image charge located at is the separation of interfaces and describes the jump in dielectrics . The function of the inverse screening length is the average between and , i.e.,

By subtracting the free-space Green’s function and taking the self Green function limit, and the self energy is then found, for

(20) |

where the first term is the contribution from the local ions, and the two series are the contribution from image charges. It can be observed that, if then and we notice that the image charges are repulsive to the ions. If then the image charges with odd are attractive, but those of even are repulsive.

The WKB formulation (20) is derived from the perturbative expansion with small , and is also accurate for large because the strong screening leads to weak image potentials. Thus the formulation is aysmptotically exact for both small and large limits. For the intermediate , the WKB approximation can be considered as an interpolation of the two limits, which, however, may produce poor prediction for the self energy of ions near the interfaces due to inaccurate estimation of the local contribution . In the following, we propose an improved approximation.

We remain the form of the image charge series, but modifying the local contribution by the similar technique of Born series (truncated at the first order) widely known in the field of quantum scattering. Consider the following two Green’s functions and which satisfy,

(21) | |||

(22) |

Let , then can be approximately,

(23) |

Subtracting Eq. (21) by Eq. (22), we have the following approximate equation for ,

(24) |

where in the right side the assumption is applied by considering is a perturbation of .

Since Eq. (24) is homogeneous in and . Now we take polar-coordinate Fourier transform in these two coordinates to get,

(25) |

where . The solution of this equation for can be expressed in an integral form by using one-dimensional Green’s function,

(26) |

Since , and is a big quantity, we have the asymptotic,

(27) |

Let , then we obtain explicit expression of the inverse Fourier transform, by setting and , which is,

(28) |

where

(29) | |||||

and

(30) |

Here, is the Bessel function, and is the incomplete Gamma function.

Then we could replace in Eq. (20) by to obtain an improved version of the WKB approximation,

(31) |

and now the image-charge effect is included. The improvement from to lies in taking into account the edge effect near the boundaries, i.e., the solvent is in confinement between two interfaces. With it, the ion interaction near interfaces is strengthened due to the weaker screening.

#### iii.3.2 Finite difference approximation

The numerical solution of the Green’s function is usually difficult due to high dimensionality – it is a function of both the source and field coordinates. We will use the properties of the geometric symmetry of the considered system and only the self Green’s function being required. These two properties allow us to develop very efficient solver to get the self energy, which has been coupled with the modified Poisson-Boltzmann equation to simulate equilibrium charged systems Xu and Maggs (2013). We briefly overview the algorithm below.

In order to reduce the dimensions, the polar symmetric Fourier transform with respect to and is first applied. Let and be the frequency correspondences of and , and let be the frequency. The Fourier transform of the generalized DH equation is a two-dimensional equation of and for each ,

(32) |

Similarly, we perform Fourier transform for the free-space Green’s function equation, , for which we get,

(33) |

To numerically solve the equations, we discretize the derivatives by central differences for and . The Dirac delta function is approximated by Kronecker delta, It should be noticed that the free-space Green’s function should be numerically approximated in order to cancel the numerical singularity of the Green’s function. Then for each , the frequency Green’s functions are obtained by finding the inverse of the coefficient matrices of linear algebra system by discretizing Eq. (32). Then the image and correlation self energy is given by the inverse Fourier transform,

(34) |

where and are diagonal elements of the approximate frequency Green’s functions. This equation is approximated by numerical integration.

## Iv Numerical results

In this section, we present numerical results by performing calculations of monovalent electrolytes. Two species of ions with valence are included, and initially the dimensionless ion densities set . The separation between the dielectric interface and the boundary takes , and thus . In all calculations, the parameters and electrode voltage (correspondingly, ) are fixed, and different strengths of and are investigated. Since the Debye length (typically, ) is the characteristic thickness of the electric double layer, this choice of does not introduce much interaction between two double layers and continuum theory is thus considered adequate to describe the charge dynamics Daiguji (2010), and we will see that introducing self-energy contribution significantly modifies the system quantities at the length scale.

The WKB and finite-difference numerical approaches for the generalized DH equation are compared. Since we have two WKB approaches, we label the results of different approaches by “WKB1”, “WKB2” and “FDM”, where WKB1 is based on expression (20) and WKB2 is the improved version Eq. (31). Also we present the results of classical PNP, denoted by “PNP” (corresponding to ), to show the difference between the PNP and the modified model. The image charge series in WKB formulations are truncated at , large enough to ensure that the truncation error is negligible. In FDM approximation, the infinite integral Eq. (34) is approximated through a cutoff at a frequency to become an integral over and a variable transformation for . We take , and 16 quadrature points in the calculations. These parameters provide high accuracy for the integration.

### iv.1 Convergence of numerical methods

In the first example, numerical schemes for the modified PNP equations with both the WKB and FDM approximations are tested. We set the parameter and the dielectric ratio , and compute the results up to time . For a water solvent at room temperature, and , the value of means the distance between electrodes is and is also physically interesting for studies of nanoscale devices. This choice of big is convenient to observe the convergence of the schemes. The dielectric constant is typical for membranes or other materials. The time step takes with varying space grid sizes and .

Fig. 2 presents the results of the spatial distribution of the net charge density, , with the three approaches and the four grid sizes. It is shown that the error of the maximum values between and is less then . By taking the results of the finest grid () as the reference, it is observed that the error is also decreased in a factor of if is halved, demonstrating all approaches are self-convergent with the second-order accuracy.

The three approaches all predict a non-monotonic curves near the interfaces. Near the cathode, the net charge is first raising to a maximum, then monotonically decays to zero at . This is due to thus the self-image charges repel the mobile ions, and in agreement with Monte Carlo simulations Messina (2002); Gan and Xu (2011); Wang and Ma (2009). In contrast, The PNP ignores the polarization effect, and thus always predicts monotonic net charge density (shown in next section). By comparing the three approaches, the charge density predicted from the WKB1 is much higher than those from the FDM and the WKB2. The maximums of the solid black lines of Fig. 2 (a)-(c) are 1.54, 1.32, and 1.30, respectively. As the FDM with the fine mesh is considered very accurate, it is concluded that the WKB1 overestimates the charge density near the interface, and lowers down the induced-charge effect. The WKB2 curves are in good agreements with the FDM results, showing a high accuracy of the analytical approximation.

### iv.2 Effect of self-energy strength

To investigate the effect of different , two groups of parameters are adopted: four self-energy strengths and 0.2, and four time snapshots and 10. When , the generalized DH equation is switched off and the model reduces to the classical PNP, and thus the PNP results are also compared to investigate how the self energy influences the results. The space grid size takes . Other parameters remain the same as the previous example.

Fig. 3 plots the net-charge-density results by the PNP, WKB2 and FDM for the four time snapshots. Again, both WKB2 and FDM are in good agreement. It is observed that image repulsion is strengthened with the increase of . For , there is an obvious maximum at . Far-field curves are overlapping at equilibrium state ().

We should see that the WKB1 approximation is already inaccurate in the case of small . To make a further comparison discussion, we plot the results of deviating from the PNP by the three approaches in Fig. 4. It is obvious the image charge should not be neglected as the charge density near the boundaries will be much smaller when it is present. One can also see that, the WKB2 prediction has been greatly improved from WKB1, and agrees well with the FDM, though the WKB2 uses the asymptotic expression for approximating the self energy.

### iv.3 Effect of dielectric ratio

Now, we investigate charge dynamics for different dielectric ratios, and We take for the modified PNP equations, and calculate the total diffusion charge in the left half of the electrolyte (near the cathode),

(35) |

as a time function. Fig. 5 presents the PNP, WKB2 and FDM results.

When , the image repulsion by FDM and WKB2 reduces the total diffusion charge, compared to the PNP. It is not difficult to understand that the image charge repulsion leads to a reduction of total charge near the surface. The FDM is slightly smaller than the WKB2 prediction. In the case of where the polarization effect is gone, the difference between the FDM (or WKB2) and the PNP is minor, and thus the electrostatic correlation between ions is weak. In the case of , the results illustrate a strong image attraction, which leads to counterion condensation on the surfaces, thus bigger total charges are predicted in comparison to the two previous cases. When time is small, the FDM and WKB2 agree well, but the deviation may be high for large and large . This demonstrates the asymptotic approach may be quantitatively inaccurate when the interfacial counterion density has a sharp change.

The FDM results also show that the total charge keeps increasing for large (or decreasing for small ) after a long time, demonstrating the existence of larger time scale in the modified PNP model. A possible explanation is the strong image attraction causes a high counterion density at surface where ion-ion correlation leads to a new relaxation time scale. For large , the solution may blow up at a certain time without including the excluded volume effect.

## V Concluding remarks

In summary, we have proposed a modified PNP model by coupling the PNP equation with the generalized DH equation and developed efficient WKB and numerical methods for the self energy. We show by numerical examples both methods are accurate, and the analytical WKB approximation is in good agreement with the difference method.

It has been investigated through simulations that image charges are involved in many many-body phenomena such as charge inversion Messina (2002); Gan et al. (2012) and like-charge attraction Wang and Wang (2013); Xu (2013). We show that the modified model can correctly predict the image-charge effect on dynamics of mobile ions. However, this work does not pay much attention on the ion-ion correlation, though this effect has been included in the modified PNP model. It is due to that the ion correlation should take into account the size effect of ions, which gives rise to a high difficulty in numerical implementation. Without accounting for the ion size effect, the solution of the modified PNP equation is unstable with the increase of applied surface voltages. The issue of overcoming this challenge is certainly our objective of future work.

## Acknowledgments

The authors acknowledge the financial support from the Natural Science Foundation of China (Grant Numbers: 11101276, and 91130012), Youth Talents Program by Chinese Organization Department, and the HPC center of SJTU. The authors thank Prof. Bob Eisenberg for careful reading and comments.

## References

- Schoch et al. (2008) R. B. Schoch, J. Han, and P. Renaud, Rev. Mod. Phys. 80, 839 (2008).
- Daiguji (2010) H. Daiguji, Chem. Soc. Rev. 39, 901 (2010).
- French et al. (2010) R. H. French, V. A. Parsegian, R. Podgornik, R. F. Rajter, A. Jagota, J. Luo, D. Asthagiri, M. K. Chaudhury, Y.-M. Chiang, S. Granick, et al., Rev. Mod. Phys. 82, 1887 (2010).
- Levin (2002) Y. Levin, Rep. Prog. Phys. 65, 1577 (2002).
- Walker et al. (2011) D. A. Walker, B. Kowalczyk, M. O. de la Cruz, and B. A. Grzybowski, Nanoscale 3, 1316 (2011).
- Messina (2009) R. Messina, J. Phys. Condens. Matter 21, 113102 (2009).
- Wan et al. (2014) L. Wan, S. Xu, M. Liao, C. Liu, and P. Sheng, Phys. Rev. X 4, 011042 (2014).
- Messina (2002) R. Messina, J. Chem. Phys. 117, 11062 (2002).
- Hatlo and Lue (2008) M. M. Hatlo and L. Lue, Soft Matter 4, 1582 (2008).
- Jho et al. (2008) Y. S. Jho, M. Kanduč, A. Naji, R. Podgornik, M. W. Kim, and P. A. Pincus, Phys. Rev. Lett. 101, 188101 (2008).
- Wang and Ma (2010) Z. Y. Wang and Y. Q. Ma, J. Phys. Chem. B 114, 13386 (2010).
- Bakhshandeh et al. (2011) A. Bakhshandeh, A. P. dos Santos, and Y. Levin, Phys. Rev. Lett. 107, 107801 (2011).
- Diehl et al. (2012) A. Diehl, A. P. dos Santos, and Y. Levin, J. Phys.: Condens. Matter 24, 284115 (2012).
- Gan et al. (2012) Z. Gan, X. Xing, and Z. Xu, J. Chem. Phys. 137, 034708 (2012).
- Xu (2013) Z. Xu, Phys. Rev. E 87, 013307 (2013).
- Wang and Wang (2013) R. Wang and Z.-G. Wang, J. Chem. Phys. 139, 124702 (2013).
- Zwanikken and de la Cruz (2013) J. W. Zwanikken and M. O. de la Cruz, Proc. Nat. Acad. Sci. USA 110, 5301 (2013).
- Chen and Eisenberg (1993) D. Chen and R. Eisenberg, Biophys. J. 64, 1405 (1993).
- Eisenberg (1996) R. S. Eisenberg, J. Membrane Biol. 150, 1 (1996).
- Eisenberg (1999) R. S. Eisenberg, J. Membrane Biol. 171, 1 (1999).
- Hollerbach et al. (2001) U. Hollerbach, D.-P. Chen, and R. Eisenberg, J. Sci. Comput. 16, 373 (2001).
- Daiguji et al. (2005) H. Daiguji, Y. Oka, and K. Shirono, Nano Letters 5, 2274 (2005).
- Markowich et al. (1990) P. Markowich, C. Ringhofer, and C. Schimeiser, Semiconductor (Springer, 1990).
- Bazant et al. (2004) M. Z. Bazant, K. Thornton, and A. Ajdari, Physical review E 70, 021506 (2004).
- Squires and Bazant (2004) T. M. Squires and M. Z. Bazant, Journal of Fluid Mechanics 509, 217 (2004).
- Singer and Norbury (2009) A. Singer and J. Norbury, SIAM J. Appl. Math. 70, 949 (2009).
- Kurnikova et al. (1999) M. G. Kurnikova, R. D. Coalson, P. Graf, and A. Nitzan, Biophys. J. 76, 642 (1999).
- Cardenas et al. (2000) A. E. Cardenas, R. D. Coalson, and M. G. Kurnikova, Biophys. J. 79, 80 (2000).
- Zheng et al. (2011) Q. Zheng, D. Chen, and G.-W. Wei, J. Comput. Phys. 230, 5239 (2011).
- Wei et al. (2012) G.-W. Wei, Q. Zheng, Z. Chen, and K. Xia, SIAM Review 54, 699 (2012).
- Flavell et al. (2014) A. Flavell, M. Machen, B. Eisenberg, J. Kabre, C. Liu, and X. Li, J. Comput. Electronics 13, 235 (2014).
- Dyrka et al. (2013) W. Dyrka, M. M. Bartuzel, and M. Kotulska, Proteins: Structure, Function, and Bioinformatics 81, 1802 (2013).
- Lockett et al. (2008) V. Lockett, R. Sedev, J. Ralston, M. Horne, and T. Rodopoulos, The J. Phys. Chem. C 112, 7486 (2008).
- Mamonov et al. (2003) A. B. Mamonov, R. D. Coalson, A. Nitzan, and M. G. Kurnikova, Biophys. J. 84, 3646 (2003).
- Busath et al. (1998) D. D. Busath, C. D. Thulin, R. W. Hendershot, L. R. Phillips, P. Maughan, C. D. Cole, N. C. Bingham, S. Morrison, L. C. Baird, R. J. Hendershot, et al., Biophys. J. 75, 2830 (1998).
- Kékicheff and Spalla (1995) P. Kékicheff and O. Spalla, Phys. Rev. Lett. 75, 1851 (1995).
- Corry et al. (2000) B. Corry, S. Kuyucak, and S.-H. Chung, Biophys. J. 78, 2364 (2000).
- Storey et al. (2008) B. D. Storey, L. R. Edwards, M. S. Kilic, and M. Z. Bazant, Physical Review E 77, 036317 (2008).
- Kilic et al. (2007) M. S. Kilic, M. Z. Bazant, and A. Ajdari, Physical Review E 75, 021503 (2007).
- Horng et al. (2012) T.-L. Horng, T.-C. Lin, C. Liu, and B. Eisenberg, J. Phys. Chem. B 116, 11422 (2012).
- Lu and Zhou (2011) B. Lu and Y. Zhou, Biophys. J. 100, 2475 (2011).
- Eisenberg et al. (2010) B. Eisenberg, Y. Hyon, and C. Liu, J. Chem. Phys. 133, 104104 (2010).
- Liu et al. (2012) W. Liu, X. Tu, and M. Zhang, J. Dynamics Diff. Equ. 24, 985 (2012).
- Nadler et al. (2003) B. Nadler, U. Hollerbach, and R. S. Eisenberg, Phys. Rev. E 68, 021905 (2003).
- Corry et al. (2003) B. Corry, S. Kuyucak, and S.-H. Chung, Biophys. J. 84, 3594 (2003).
- Graf et al. (2004) P. Graf, M. G. Kurnikova, R. D. Coalson, and A. Nitzan, J. Phys. Chem. B 108, 2006 (2004).
- Hyon et al. (2011) Y. Hyon, B. Eisenberg, and C. Liu, Commun. Math. Sci 9, 459 (2011).
- Luo et al. (2006) G. Luo, S. Malkova, J. Yoon, D. G. Schultz, B. Lin, M. Meron, I. Benjamin, P. VanÃ½sek, and M. L. Schlossman, Science 311, 216 (2006).
- Podgornik (1989) R. Podgornik, J. Chem. Phys. 91, 5840 (1989).
- Netz and Orland (2000) R. R. Netz and H. Orland, Eur. Phys. J. E 1, 203 (2000).
- Netz and Orland (2003) R. R. Netz and H. Orland, Eur. Phys. J. E 11, 301 (2003).
- Bonthuis et al. (2011) D. J. Bonthuis, S. Gekle, and R. R. Netz, Phys. Rev. Lett. 107, 166102 (2011).
- Bonthuis and Netz (2013) D. J. Bonthuis and R. R. Netz, The J. Phys. Chem. B 117, 11397 (2013).
- Wang (2010) Z. G. Wang, Phys. Rev. E 81, 021501 (2010).
- Cherepanov et al. (2003) D. A. Cherepanov, B. A. Feniouk, W. Junge, and A. Y. Mulkidjanian, Biophys. J. 85, 1307 (2003).
- Avdeev and Martynov (1986) S. M. Avdeev and G. A. Martynov, Colloid J. USSR 48, 535 (1986).
- Buyukdagli et al. (2012) S. Buyukdagli, C. V. Achim, and T. Ala-Nissila, J. Chem. Phys. 137, 104902 (2012).
- Yaroshchuk (2000) A. E. Yaroshchuk, Adv.Colloid Interface Sci. 85, 193 (2000).
- Lau et al. (2002) A. W. C. Lau, D. B. Lukatsky, P. Pincus, and S. A. Safran, Phys. Rev. E 65, 051502 (2002).
- Buff and Stillinger (1963) F. P. Buff and F. H. Stillinger, J. Chem. Phys. 39, 1911 (1963).
- Xu and Maggs (2013) Z. Xu and A. C. Maggs, J. Comput. Phys. (2014), eprint http://dx.doi.org/10.1016/j.jcp.2014.07.004
- Lin et al. (2011a) L. Lin, C. Yang, J. C. Meza, J. Lu, L. Ying, and W. E, ACM Trans. Math. Softw. 37, 40:1 (2011a).
- Lin et al. (2011b) L. Lin, C. Yang, J. Lu, L. Ying, and W. E, SIAM J. Sci. Comput. 33, 1329 (2011b), eprint http://epubs.siam.org/doi/pdf/10.1137/09077432X.
- Pasquali and Maggs (2009) S. Pasquali and A. C. Maggs, Phys. Rev. A 79, 020102 (2009).
- George (1973) A. George, SIAM J. Numer. Anal. 10, 345 (1973).
- Davis (2006) T. A. Davis, Direct Methods for Sparse Linear Systems (SIAM, Philadelphia, 2006).
- Gan and Xu (2011) Z. Gan and Z. Xu, Phys. Rev. E 84, 016705 (2011).
- Wang and Ma (2009) Z. Y. Wang and Y. Q. Ma, J. Chem. Phys. 131, 244715 (2009).