{'fraction_non_alphanumeric': 0.03889914815035444, 'fraction_numerical': 0.013284088878298683, 'mean_word_length': 4.546836279530811, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Self-energy-modified Poisson-Nernst-Planck equations: WKB approximation and finite-difference approaches

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\"uckel (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.


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 [1, 2]; 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 ℓ B , and the Debye length of the electrolyte ℓ D ), complex electrostatic phenomena such as the charge inversion and the likecharge attraction have been often observed experimentally. These phenomena are beyond mean-field theoretical explanations and motivate many challenges to computational and modelling communities [3].
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 [3][4][5][6][7]. 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 correla-tion 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 [8][9][10][11][12][13][14][15][16][17] 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 [18][19][20][21]. The system of equations has also been widely applied in continuum theories of narrow channels by coupling with the Stokes equations to model nanofluidics [2,22]. The PNP equations are also known as "drift-diffusion equations" and have much use in semiconductors [23]. Asymptotic analysis on the PNP equations has been made to understand properties of charge diffusion of electrochemical energy systems [24][25][26]. From the aspect of numerical analysis, various numerical methods have been proposed [27][28][29][30][31][32] 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 [33]. Analysis of current-concentration curves [34] also showed significantly different trends between the experimental data [35] 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 [36]. The comparison between the PNP theory and particle-based Brownian dynamics simulations [37] for cylindrical pores of varying sizes demonstrated the invalidity of the meanfield theories when the cylindrical radius is less than 2ℓ D , 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 [38][39][40][41][42][43] or the dielectric self energy [44][45][46][47]. 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., [16,48]. 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 [49][50][51] 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, ǫ = ℓ D /L and q = ℓ B /L, where L 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 L 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 [16,17]. 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 [52,53]. 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 U i of the ions, where c i is the concentration of the ions of species i, D i is its diffusion constant, β = 1/k B T where k B is the Boltzmann constant and T is the temperature.
In the mean-field description of electrostatics, the potential energy U i takes the mean potential energy U i = z i eΦ, and the electric potential Φ is governed by the Poisson equation, where z i is the valence, ρ f e is the fixed charge, ε 0 is the vacuum dielectric permittivity and ε(r) 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 inducedcharge 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 U i can be expressed as the mean potential energy plus a correction, where u i is the self energy of a unit test ion of the ith transported species. Given the mean potential Φ, an accurate description of u i 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 [48], or by a reaction-field formulation [45]. In the self-consistent Gaussian field approximation [50,51,54], this quantity can be defined through a Green's function G(r, r ′ ), described by the generalized Debye-Hückel (DH) equation, and u i is then the self Green's function limit, is the Green's function in free space to remove the invariable singularity, and I 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 ε ′ i describes the function has been locally modified due to the excluded volume of the test ion, where a i is the ionic radius, and ε eff 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), I is expressed as, This modification is particularly useful to remove the self-energy singularity when the permittivity of solvent is space-dependent [53,55], but greatly increases numerical difficulty. To avoid this difficulty, we will not discuss this size effect in this work, and simply take the limit a i → 0, 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 c i0 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, Together with the generalized DH equation (4), it has been studied by Avdeev and Martynov [56] using the Debye closure of the BBGKY chain, Netz and Orland [50,51] using Gaussian variational field theory, and in much recent work (see [9,[57][58][59] 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 x = ±L (see Fig. 1), which is a simple model of electrochemical systems. We study the case of sharp dielectric permittivity ε which takes the water permittivity ε W for |x| < (1 + ξ)L with a small ξ > 0, and the alternative value ε B 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.
FIG. 1: (Color online) Schematic illustration of the electrolyte between two dielectric interfaces with applied voltages. When εB < εW , the coefficient γ > 0 and the mobile ions near the interfaces are repelled by image charges.

A. Dimensionless equations, boundary and initial conditions
Let ℓ B = βe 2 /(4πε 0 ε W ) be the Bjerrum length at the water solvent, at which distance the interaction energy of two unit point charges is k B T . We take L as the reference length scale, D 0 as a typical diffusion constant, and c 0 as a typical ion condensation. We assume ions have uniform diffusion constant and define the another length scale ℓ D = 1/ √ 8πℓ B c 0 . We note, when c 0 is the bulk ion concentration for symmetric monovalent salt, ℓ D is the Debye screening length.
The following dimensionless parameters and variables are defined: We drop the tildes of all new variables. Then the modified PNP and generalized DH equations are, where ǫ = ℓ D /L and q = ℓ B /L are two length-scale parameters. We see ǫ and 2ǫ 2 describe the effective diffusion constant of the Nernst-Planck equation, and dielectric permittivity of the Poisson and generalized DH equations, and q 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 −1 < x < 1, while the generalized DH equation is defined on the whole threedimensional domain with implied interface conditions at x = ±(1 + ξ). We assume completely blocking electrodes [24] imposed with a time-varying external potential, V ± in dimensionless unit, on the electrodes. Then, at x = ±1, 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, In this work, we take η = 0 to use the Dirichlet condition for the potential. The Robin boundary condition, i.e., η = 0, is often used to account for the dielectric electrodes in literature [24,31] by considering the tight counterion adsorption at the Stern layer. In the Green's function equation, the dielectric permittivity is discontinuous at interfaces at x = ±(1 + ξ). We find as a result of using a positive ξ, the higher-order singularity of the self energy can be removed when x ′ → ±1. The dimensionless dielectric permittivity is then ε = 1 for |x| ≤ (1 + ξ), and ε B /ε W elsewhere. Initially, it is assumed that the ions are in equilibrium without applying the external potential which starts to function at t = 0. Therefore, we can set constant initial ionic concentrations, c i (x, 0) = c 0 i , under electricneutral constraint i z i c 0 i = 0. For monovalent salt we are studying, c 0 i = 1 for i = 1 and 2. In addition, we assume there is no fixed charge ρ f = 0 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 q.

B. Discretization of the PNP equations
The Nernst-Planck and the Poisson equations are approximated by finite-difference discretizations. Let ∆t and ∆x be the time and space steps, and denote the nth time and kth space grids by t n = n∆t and x k = −1 + k∆x, respectively. The second-order semi-implicit time-stepping scheme is employed for the Nernst-Planck equation, where U i = z i Φ + 1 2 z 2 i u, and c Here the scheme of diffusion term is a central difference, and that of the advection term is an extrapolation for approximating U n+ 1 2 i . By this discretization, we gain the secondorder 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 ∆t ∝ ∆x. Since it is a three-point scheme in time, we should use the backward scheme for U i at the initial step n = 0.
Concerning the spatial discretization, because the dielectric function ε = 1 is constant in interval −1 < x < 1 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 c i and U i 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 c n+1 i , the mean potential at (n + 1)th time step can be simply obtained by central discretization with the boundary condition.

C. 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, where the generalized inverse Debye length, κ(x) = i z 2 i c n+1 i /(2ǫ 2 ) between two interfaces and zero otherwise. Here we set q = 1 without loss of generality. Solving this equation gives the Green's function at time t n+1 . 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 [16,57], which avoids numerical solution of the high-dimensional problem. We propose a simple WKB expression by improving the idea of Buff and Stillinger [60] which has been used recently for studying the double layer interaction by Wang and Wang [16]. 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.

WKB approximation
In the WKB approximation introduced by Buff and Stillinger [60] for one-interface problem, the Green's function is first exactly found in the case of κ(x) 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, where r ℓ is the distance between r and the ℓth image charge located at ℓD + (−1) ℓ x ′ , y ′ , z ′ , D = 2(1 + ξ) is the separation of interfaces and γ describes the jump in dielectrics γ = εW −εB εW +εB . The function of the inverse screening length is the average between x and x ′ , i.e., By subtracting the free-space Green's function and taking the self Green function limit,κ(x, x ′ ) = κ(x) and the self energy u is then found, for −1 ≤ x ≤ 1, where the first term u loc = −κ(x) is the contribution from the local ions, and the two series are the contribution from image charges. It can be observed that, if ε B < ε W , then γ > 0 and we notice that the image charges are repulsive to the ions. If ε B > ε W , 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 −u loc . In the following, we propose an improved approximation.
We remain the form of the image charge series, but modifying the local contribution κ(x) = − u loc 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 G 1 (r, r ′ ) and G 0 (r, r ′ ) which satisfy, Let δG = G 1 − G 0 , thenũ loc can be approximately, Subtracting Eq. (21) by Eq. (22), we have the following approximate equation for δG, where in the right side the assumption G 1 = e −κ(x)|r−r ′ | /|r − r ′ | is applied by considering G 1 is a perturbation of G 1 . Since Eq. (24) is homogeneous in y and z. Now we take polar-coordinate Fourier transform in these two coordinates to get, where b(x) = ω 2 + κ(x) 2 . The solution of this equation for x ∈ [−D/2, D/2] can be expressed in an integral form by using one-dimensional Green's function, Since ǫ ≪ 1, κ ≫ 1 and b is a big quantity, we have the asymptotic, Let ρ = (y − y ′ ) 2 + (z − z ′ ) 2 , then we obtain explicit expression of the inverse Fourier transform, by setting y ′ → y and z ′ → z, which is, where and Here, J 0 is the Bessel function, and Γ[α, z] . = ∞ z t α−1 e −t dt is the incomplete Gamma function. Then we could replace κ(x) in Eq. (20) by κ(x) = − u loc to obtain an improved version of the WKB approximation, 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.

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 [61]. We briefly overview the algorithm below.
In order to reduce the dimensions, the polar symmetric Fourier transform with respect to y − y ′ and z − z ′ is first applied. Let G and G 0 be the frequency correspondences of G and G 0 , and let ω be the frequency. The Fourier transform of the generalized DH equation is a two-dimensional equation of x and x ′ for each ω, (32) Similarly, we perform Fourier transform for the free-space Green's function equation, −∇ 2 G 0 = 4πqδ(r − r ′ ), for which we get, To numerically solve the equations, we discretize the derivatives by central differences for G and G 0 . The Dirac delta function is approximated by Kronecker delta, δ(x j − x k ) ≈ δ jk /∆x. It should be noticed that the freespace 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, where G jj (ω) and G 0,jj (ω) are diagonal elements of the approximate frequency Green's functions. This equation is approximated by numerical integration. We remark that the direct inversion of a matrix is expensive, and we should apply some sparse inversion technique to find the diagonals of an inverse matrix; for literature, see [62][63][64][65][66].

IV. NUMERICAL RESULTS
In this section, we present numerical results by performing calculations of monovalent electrolytes. Two species of ions with valence z 1,2 = ±1 are included, and initially the dimensionless ion densities set c 1,2 = 1. The separation between the dielectric interface and the boundary takes ξ = 0.06, and thus D = 2.12. In all calculations, the parameters ǫ = 0.2 and electrode voltage V ± = ±1 (correspondingly, ±25.6mV ) are fixed, and different strengths of q and ε B /ε W are investigated. Since the Debye length ℓ D (typically, ∼ 1 − 100nm) 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 [2], 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 q = 0), to show the difference between the PNP and the modified model. The image charge series in WKB formulations are truncated at ℓ = 10, 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 [0, Ω] and a variable transformation ω = e µv − 1 for µ > 0. We take µ = 1, Ω = 1024 and 16 quadrature points in the calculations. These parameters provide high accuracy for the integration.

A. 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 q = 0.2 and the dielectric ratio ε B /ε W = 0.05, and compute the results up to time T = 2. For a water solvent at room temperature, ε W ≈ 80 and ℓ B ≈ 0.7nm, the value of q = 0.2 means the distance between electrodes is L ∼ 3.5nm and is also physically interesting for studies of nanoscale devices. This choice of big q is convenient to observe the convergence of the schemes. The dielectric constant ε B = 4 is typical for membranes or other materials. The time step takes ∆t = ∆x/2 with varying space grid sizes ∆x = 1/100, 1/200, 1/400 and 1/800. Fig. 2 presents the results of the spatial distribution of the net charge density, z 1 c 1 (x, T ) + z 2 c 2 (x, T ), with the three approaches and the four grid sizes. It is shown that the error of the maximum values between ∆x ≤ 1/400 and 1/800 is less then 2%. By taking the results of the finest grid (∆x = 1/800) as the reference, it is observed that the error is also decreased in a factor of ∼ 4 if ∆x 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 x = 0. This is due to ε B < ε W thus the self-image charges repel the mobile ions, and in agreement with Monte Carlo simulations [8,67,68]. 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.

B. Effect of self-energy strength
To investigate the effect of different q, two groups of parameters are adopted: four self-energy strengths q = 0, 0.05, 0.1 and 0.2, and four time snapshots T = 0.2, 0.5, 2 and 10. When q = 0, 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 1/800. 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 q. For q = 0.2, there is an obvious maximum at x ≈ −0.95. Far-field curves are overlapping at equilibrium state (T = 10).
We should see that the WKB1 approximation is already inaccurate in the case of small q. 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.

C. Effect of dielectric ratio
Now, we investigate charge dynamics for different dielectric ratios, ε B /ε W = 1/20, 1 and 20. We take q = 0.1 for the modified PNP equations, and calculate the total diffusion charge in the left half of the electrolyte (near the cathode), as a time function. Fig. 5 presents the PNP, WKB2 and FDM results. When ε B /ε W < 1, 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 ε B /ε W = 1 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 ε B /ε W > 1, 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 T and large ε B . 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 ε B (or decreasing for small ε B ) 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 ε B , 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 [8,14] and like-charge attraction [15,16]. 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.