Energy variational approach to study charge inversion (layering) near charged walls

We introduce a mathematical model, which describes the charge inversion phenomena 
in systems with a charged wall or boundary. 
This model may prove helpful in understanding semiconductor devices, ion channels, 
and electrochemical systems like batteries that depend on complex distributions of charge for their function. 
The mathematical model is derived using the energy variational approach that 
takes into account ion diffusion, electrostatics, finite size effects, and 
specific boundary behavior. 
In ion dynamic theory, a well-known system of equations is the Poisson-Nernst-Planck (PNP) equation 
that includes entropic and electrostatic energy. 
The PNP type of equation can also be derived 
by the energy variational approach. 
However, the PNP equations have not produced the charge inversion/layering 
in charged wall situations presumably 
because the conventional PNP does not include the finite size of ions 
and other physical features needed to create the charge inversion. 
In this paper, we investigate the key features 
needed to produce the charge inversion phenomena using a mathematical model, 
the energy variational approach. 
One of the key features is a finite size (finite volume) effect, 
which is an unavoidable property of ions important for their dynamics on small scales. 
The other is an interfacial constraint to capture the spatial variation of electroneutrality 
in systems with charged walls. 
The interfacial constraint is established 
by the diffusive interface approach that approximately describes the boundary effect 
produced by the charged wall. 
The energy variational approach gives us a mathematically self-consistent way to introduce the interfacial constraint. 
We mainly discuss those two key features in this paper. 
Employing the energy variational approach, we derive a non-local partial differential equation 
with a total energy consisting of the entropic energy, electrostatic energy, 
repulsion energy representing the excluded volume effect, 
and the contribution of an interfacial constraint related to overall electroneutrality 
between bulk/bath and wall. 
The resulting mathematical model produces the charge inversion phenomena near charged walls. 
We compare the computational results of the mathematical model to those of Monte-Carlo computations.


1.
Introduction. Ionic solutions are in many ways the most important mixtures people, and thus mathematicians, encounter. All of life occurs in ionic solutions; much of chemistry is done in ionic solutions, all our thoughts are nerve cells working in ionic solutions. Understanding ionic solutions has been entwined with understanding life since ions were discovered and their fluxes investigated by Fick, a biologist, actually a physiologist [35,38,51]. We deal with a particular phenomenon that arises from the interaction of ions and a charged wall. Ions have charge and finite diameter. They do not overlap. Thus, when ions are near a charged wall, complex distributions occur if the ions have different charge and/or diameter. These complex distributions near the walls have historically been important in characterizing non-ideal behavior of solutions and theories. They also have technological importance as the source of the strength of cement [24,41,58]. The complex distributions are closely related to distributions in and near proteins and lipid membranes [16,18,20].
Complex distributions of charge are important in many applications. The switching and amplifying characteristics of transistors, and thus integrated circuits, arise in depletion layers of charge [46,52,62,64]. The selectivity properties of cardiac calcium channels [6,9], voltage activated sodium channels of nerve [8,22], and Ryanodine Receptors [27,28,29,30,32] arise from distributions of mobile and permanent charge inside proteins [19,47,48]. Batteries are just one example of electrochemical devices that depend on complex distributions of charge near electrodes. Complex distributions of charge are often layers, sometimes layers with alternating signs of net charge. The resulting distributions of potential are of course even more complex and can produce current voltage relations with maxima and minima, and time dependent behavior of active 'negative conductance' devices.
Charge inversion (layering) is an alternating accumulation of matter and charge near a charged wall, surface, or boundary. When the charged wall is exposed to electrolyte solution, the electric force (and sometimes potential) reverses polarity in some regions due to the complex spatial dependence of the accumulation of counterions near the wall [4]. We introduce a mathematical model to capture the essential features of charge inversion (layering) phenomena in the presence of a charged wall/boundary. The mathematical model is established by the energy variational approach [23,36,37,59]. The energy variational approach has been previously used to compute transport in homogeneous systems. The resulting system of partial differential equations is the Poisson-Nernst-Planck (PNP) equation, which has well-known applications in semiconductor problems [40,45,46,49,60]. Historical references to the chemical and biophysical literature on PNP can be found in [3,17,21]. The classical PNP equation does not include the excluded volume effects produced by the finite diameter of ions; It can be proven that it does not produce layering behavior for that reason. The conventional PNP equation, and also classical Poisson Boltzmann (PB) equation for that matter, are only valid for diluted solutions. In fact, under most physical boundary conditions, it is not possible to satisfy such diluted conditions everywhere in for bulk diluted solutions, since ions can (and will) accumulate near the boundaries, as mentioned before in situations like batteries and ion channels. In [23,37], we introduced a modified PNP including finite diameter of ions. We want to point out that the specific way of including the finite size effects is only the first step to address the physical nondilute situations of ion particles. However, all the physics (and chemistry) can be included in the nonlocal energy functional. Once the functional is determined, one can uniquely derive the resulting system using the energetic variational approach [23]. From the point of view analysis and numerics, the methods to deal with the nonlocal terms in the resulting partial different equations are also crucial for any more realistic treatment of size effects. Finally, we notice that some 1-dimensional simulations do not deal robustly with layering found in simulations in high dimensional spaces (in 2, 3-dimensional spaces), although it can reproduce some parts of the phenomena [56,57]. For that reason, we seek a key feature/factor that would allow a modified PNP to recover the charge inversion phenomena robustly.
The key idea for a new PNP equation comes from the observation of charge density profiles in the double layer region near the charged wall in a large number of computations using both Monte Carlo (MC) simulations and the PNP system. It seems (qualitatively) that charge neutrality is enforced too strongly in the equation of [23,37] in one dimension. It seemed likely to us that the condition of bulk electroneutrality needs to be relaxed close to the wall to get a robust one dimensional treatment of layering, comparable to simulation results in one and higher dimensions.
Specifically, comparisons between modified PNP [23,37] and MC results of spheres near a charged wall were done in a wide variety of conditions. When the wall had small charge, the numerical results of PNP with excluded volume and MC show a good agreement when solutions contained monovalent ions or mixtures of monovalent and divalent ions. But large wall charge situations show quite different behaviors near the charged wall. MC results show a charge inversion both of charge density and electrostatic potential, but modified PNP does not show a charge inversion of the electrostatic potential. For that reason, we introduce a new modified PNP that we call PNP-FS-IF because it includes (i ) a nonlocal contribution of finite size (FS) and (ii ) an interfacial constraint (IF) related to an electroneutrality condition near charged walls/boundaries.
We relax the constraint of bulk electroneutrality near the charged wall using the diffusive interface approach. The diffusive interface approach is a well-known mathematical procedure to describe interactions between different phases and dynamics of interfaces [1,10,11,12,43,67] that apparently has not been used before to deal with charge layering. The diffusive interface approach allows us to establish an interfacial constraint for the relaxation of bulk electroneutrality as it approaches the wall, and effectively gives us a description of the influence of the wall charge on the ionic solution. However, establishing a system containing the interfacial constraint is another task. Since electroneutrality is a principal principle of ionic solutions, and an important but delicate feature of their physics, the mathematics must be done with great care to be consistent. More specifically, tiny changes in net charge anywhere in a solution have a large effect on the potential 'everywhere'. The very small deviation from equality in concentration of negative and positive charge is exactly the change in charge that controls the electrical forces. The electric forces are large -often much larger than the effective force of diffusion-even though they arise from such tiny fractions of all the ions present (typically, 10 −15 ). An inconsistent procedure to put electroneutrality constraint into a system of equations could easily break the electroneutrality of the system.
The energy variational approach gives us a self-consistent way to apply an interfacial constraint related to electroneutrality. We add a penalization 'energy' using the interfacial constraint to define a new total energy. Then, we apply the energy variational approach to this new energy and derive a new modified equation, PNP-FS-IF. In this paper, we shall show a numerical verification of the PNP-FS-IF equations only in one dimensional space. The PNP-FS-IF formulation may prove helpful in understanding semiconductor devices, ion channels, and electrochemical systems like batteries that depend on complex distributions of charge for their function.
The organization of the paper is the following: First, we discuss MC computations related to the charge inversion phenomena. In the next section, we present the derivation of a modified PNP equation including the contribution of the excluded volume effect and the interfacial constraint, using the the energetic variational approach. In section 4, we present several numerical experiments as the verification of the resulting system of partial differential equations. Finally we propose PNP-FS-IF as a useful model of interfacial charge accumulation in a variety of applications, from electrodes and batteries to lipid bilayers and ion channels.

2.
Monte-Carlo methods. The MC computations provide benchmarks for the charge inversion phenomena [31]. Even though the MC computation is not identical to our other numerical methods, in order to show the capability of PNP-FS-IF model to capture the charge inversion in phenomenal viewpoints, we perform the MC simulations to compare with the energetic variational approach results. In the simulation model, ions are represented as spheres with Pauling radii and a point charge at their centers. The system is contained on two sides by parallel walls, which may have a continuous charge, given in C/m 2 . Hard interactions mean that no overlap is allowed between ions or between an ion and a wall. Water is described in the "primitive model" as implicit solvent with a continuum dielectric with a dielectric coefficient of 78.5. The temperature of the system was 298.15K.
Coulombic interactions are summed for every ion-pair to determine the energy of the system. For ions of species i and j, the interaction potential is where r is the distance between the ion centers, z i is the valence of ith ion, q is the fundamental charge, and ε is the permittivity of free space. a ij is the average radius of ions of species i and j.
The walls are represented as infinite sheets of charge. The electrostatics are modeled with Torre-Valleau boundary conditions [63]. The interaction between an ion and the wall is where d i is the distance from ith ion to the closest point on the wall and σ is the surface charge density of the wall.
3. Modified Poisson-Nernst-Planck equations: Energy variational approach. In this section, we derive a mathematical model, a modified PNP equation, using the energy variational approach. We consider the entropic energy for Brownian motion of spheres, the electrostatic energy for Coulomb interactions for charged ions, the repulsive potential energy for excluded volume effects, and an interfacial constraint related to electroneutrality.
First, we present the derivation of the PNP equation including only the basic ion transport with diffusion and drift from Brownian motion and the electrostatic potential, respectively. To derive the PNP equation, we first define the total energy with the entropic and electrostatic energies, where k B is the Boltzmann constant, T is the absolute temperature, N is the number of ions, ε is the dielectric constant and c i is the charge density for ith ion species, and φ is the electrostatic potential. Since the electrostatic potential φ satisfies the Maxwell equation given by the following Poisson equation: where z i is the valence of ith ion species, ρ 0 is the permanent charge density, and q is the unit charge, the total energy E total P N P can be rewritten as with suitable boundary conditions. Now, we take a variational derivative of E total P N P with respect to the charge density c i using the explicit form of the electrostatic potential, which is the solution of (4), where G( x, y) is the Gaussian kernel. Then we have that Hence, the chemical potential µ P N P i for ith ion species is obtained by Since the ionic flux J P N P i is given by where D i is the diffusion constant for ith ion species, the charge conservation leads us to the desired equation, for i = 1, · · · , N .

YUNKYONG HYON, JAMES E. FONSECA, BOB EISENBERG AND CHUN LIU
Combining with the Poisson equation (4), we have the PNP equation The system of partial differential equations (11), (12) naturally satisfies the following energy dissipation law: This is one advantage of the energy variational approach.
In the next subsections we present/discuss the derivation of the finite size effects and then the interfacial constraint related to the electroneutrality condition between bulk and charged wall. Since the detailed derivation of finite size effect can also be found in [37], we only summarize its derivation in the next section so this paper is reasonably complete.
3.1. Excluded volume effects. The excluded volume effect is a crucial feature of ion interaction on small scales, with important effects extending into macroscopic scales, and is based on a repulsive interaction between ions, or ion and wall. An excluded volume effect is established by a mathematical model with a regularized repulsive interaction potential in [23,36]. The interaction potential is modeled by the repulsive part of the Lennard-Jones (LJ) potential. Here, we briefly present the derivation.
The energy of repulsion is defined by where the repulsive interaction potential is for ith and jth ions located at x and y with the radii a i , a j , respectively, where ε i,j is an energy constant which is symmetric and empirically chosen to match the experimental data. The variational derivative of the repulsion energy E repulsion LJ is δE repulsion Thus, the chemical potential µ repulsion i is given by We then have a new equation for charged ion transport combining the excluded volume contribution with the Poisson equation, The system (18), (19) satisfies the following energy law: 3.2. The interfacial constraint between bulk/bath and charged wall. We introduce an interfacial constraint between bulk/bath and charged wall related to electroneutrality trying to catch the essence of the nonlinear behavior of charged ionic solutions near the wall, membrane, or boundary. The bulk/bath is the region far from the charged wall. The bulk/bath is the region in Figure 1 that excludes the double layer. The electroneutrality condition is an important feature in electromagnetic theory [39] and physical chemistry [2,13,14,15,25,26,33,34,42,50,53,54,55,61]. It is a key element in our approach to establish a mathematical model for the charge inversion based on PNP theory. We employ the diffusive interface method and the energy variational approach to establish the system of equations with the interfacial constraint. We here emphasize that the diffusive interface method allows us make an effective boundary model of actions of the wall charge. The energy variational approach gives us a mathematically self-consistent way to obtain, derive, and apply an interfacial constraint for electroneutrality.
In Electrodynamics [39], when a charged wall/membrane is present, charge density and the electrostatic potential show a nonlinear and complicated behavior (double layer) near the charged wall. If the wall is highly charged, complicated phenomena arise, for instance, charge inversion (layering) which is an alternating electric force (or even potential) or alternating accumulation of charge caused by the induced electric field.
The electroneutrality condition is usually given by a pointwise constraint in addition to a global constraint on the entire system. According to double layer theory, the pointwise (local) electroneutrality condition is valid for the whole domain when no permanent charge is present. When the wall has permanent charge, the pointwise electroneutrality condition is valid in the bulk but it is not valid near the charged wall. Charge accumulation is called 'double layer capacitance' in the old chemistry literature. This is also most of the capacitance of field effect transistors and so is a dominant determinant of the speed and power of much of our digital technology both in cpu and memory chips. The pointwise electroneutrality constraint guarantees the neutrality at every point, while the global constraint guarantees neutrality for the total amount of charge in the entire system. We suspect that the pointwise electroneutrality condition may be too strong to be applied just against the wall in the charged wall situation, at least in our one dimensional calculations, and perhaps in general. Moreover, it is hard to determine the spread of the effect of wall charge into the bulk solution while enforcing pointwise electroneutrality. Relaxing point wise electroneutrality allows storage of charge and the resulting effects of charge, namely capacitance, displacement current, electric field, and electric forces. The global electroneutrality condition can be easily combined with the energy variational approach to get a mathematically self-consistent contribution to the electroneutrality condition. We consider the global constraint in the modeling of the charge layering (inversion). Then we consider a simple case in which the wall is located at the boundary of the domain Ω.
where c 0 is a nonzero constant. Since the permanent charge is located only at the boundary (wall), the electroneutrality condition for the bulk should be However, the double layer phenomena (under the influence of the permanent charge) shows that the condition (23) is not valid for the interfacial (double layer) region near the charged wall because of the influence of the electric field caused by the permanent charge. The total charge of the solution is not zero. Rather it equals the charge on the boundary. The solution as a whole has net charge. This charge is stored in what an electrically oriented physicist or engineer calls the 'capacitance to ground', or a physical chemist might call the 'self-energy'. Hence, the crucial procedure in the modeling is how to model the influence of permanent charge in the interfacial region between the boundary and the bulk. When applying the bulk electroneutrality condition (23) to the double layer theory, we suspect that the bulk electroneutrality condition should be relaxed from bulk to wall. One good mathematical tool for interfacial behaviors or interactions is the diffusive interface approach that was developed to describe the behaviors of immiscible fluids, mostly related to the Ginzburg-Landau mixing energy [1,10,11,12,43,67]. To do this, we define a function ψ representing the diffusive interface corresponding to the interfacial region between bulk and wall. Using ψ we define an interfacial constraint as Remark 2. We try to catch the essence of the nonlinear interfacial behavior near a charged wall with the phase function ψ that can be defined as a smooth function satisfying ψ = 1, in bulk, 0, at the charged wall.
The function ψ is described by the hyperbolic tangent tanh(·) as shown in Figure  1. Now we define the total energy including the interfacial constraint as a penalization, where λ(> 0) is a Lagrange multiplier that can be interpreted as a regulator or measure of the relative strength of interfacial electroneutrality in the total energy. It depends, of course, on the permanent charge ρ 0 , of the interfacial constraint on the energy. Then we take the variational derivative ofẼ total P N P with respect to c i to obtain its chemical potential, Remark 3. The determination of the thickness of the diffusive interface region for the function ψ and the Lagrange multiplier λ are a very important part of our procedure. The thickness is the width of the region where wall charge makes its influence felt, and it is directly connected to a certain region where the electrostatic potential is steeply changing near the charged wall. Hence, the thickness can be dynamically determined by examining the electrostatic potential profile obtained in PNP computations in time and comparing it to some known or desired behavior. However, the value of thickness might be fixed as an input parameter by the stationary data obtained in PNP or MC computations. These are in effect part of a model of the charge storage, capacitance, displacement current, and electric field of the double layer and charged wall.
The coefficient λ determines the strength of the interfacial constraint depending on the permanent charge ρ 0 . The value of λ can be approximately determined, through a trial and error manner, by comparison with the MC stationary data. Eventually, the methods of inverse problems may be used to determine the 'optimal' estimates and properties of the function ψ and the coefficient λ. The value of the coefficient must be large enough to produce charge inversion (layering), and must be zero when the wall is uncharged.
The system of partial differential equations for basic charge transport and the interfacial constraint is Remark 4. If the phase function ψ and λ are constants in the whole domain, then the system (27), (28) recovers the PNP equation (11), (12).
The system (27), (28) also satisfies the following energy dissipation law: Finally, combining the PNP equation with the excluded volume effect, and the interfacial constraint, we have the desired PNP-FS-IF equation as follows: Also, we easily see that the system (30) satisfies the energy dissipation law: In the next section we perform several numerical experiments for a verification of the new system of equations, PNP-FS-IF, and compare its numerical results to those of MC computations. 4. Numerical simulations. We consider two different situations. One is an uncharged wall/boundary situation, which can give a qualitative comparison between PNP-FS-IF and MC for the basic dynamics of ions. The other is a charged wall situation, in which we can observe the phenomena (charge inversion) that is the target of this paper.
We compare the charge inversion results between PNP-FS-IF including the interfacial constraint and MC computations. The computational domains are chosen as [−5, 5](nm), [−2, 2](nm) for charged and uncharged wall situations, respectively, which can be easily extended to a three dimensional domain with periodicity and free boundary arguments. These domains provide a good description of the ion concentrations and potential and their dependence on the charge of the wall, and other parameters. An inadequate bulk region may cause undesirable interaction between the walls. MC computations use an extended domain via periodic boundary conditions and take cross-sectional results. The setting of a computational domain allows us to compare those results. However, we may have a partial effect of finite size in using PNP-FS-IF.

Remark 5.
We do not know if this electroneutrality constraint would be needed in a three dimensional calculation, because three dimensional variational calculations have not yet been computed. They pose significant computational problems of efficiency and accuracy. An investigation of the excluded volume effect and the electroneutral constraint in high dimensional space will appear elsewhere. It is important to remember that low dimensional approximations will still have important advantages after three dimensional calculations are done. Low dimensional approximations will be much faster, will allow physical insight, and can be included in inverse problem calculations much more easily than the necessarily larger, slower, and clumsier three dimensional calculations. A one dimensional model with an interfacial electroneutrality constraint might prove to be a (more or less) universal model of interfacial charge accumulation, capacitance, displacement current, and electric field applicable to double layer capacitance, surface charge of lipid bilayers, and even the capacitance of field effect transistors.
For numerical computations of PNP-FS-IF, we use the edge average finite element (EAFE) method (which was developed for solving drift-diffusion equations [66]) with Dirichlet-Neumann boundary condition throughout the numerical computations. We also employ the semi-implicit time discretization to treat the nonlocal repulsion term. A computational algorithm dealing with the nonlocal repulsion for modified PNP system is found in [37]. We present a similar algorithm in Algorithm 4.1 to solve PNP-FS-IF system for self-containment. Moreover, since the kernel of repulsive potential is intrinsically singular, we compute the contribution of the nonlocal repulsive term by cutting off the domain near a singularity for numerical stability. For instance, a cut-off can be chosen as | x − y| ≥ R i,j where 0 < R i,j ≤ a i + a j for i, j = 1, · · · , N . More detailed discussion can be found in [37].
The MC simulations employ the Metropolis algorithm which allows us to sample the phase space with (thermodynamic) energies proportional to the Boltzmann factor, thereby ensuring that the simulation outputs represent the system at thermodynamic equilibrium. The thermodynamic energy is not the energy of the variational procedure.
There are four MC events: 1. moving an ion a small distance (less than its radius), 2. moving a ion to a random location within the system, 3. creating a neutral group of ions (a salt), and 4. destroying a salt. The event frequency can be tuned for computational purposes and methods of sampling have been described previously [5,6,7,44].
All MC trials are subject to an acceptance criteria. The probability, P , of accepting a MC trial allows the system to move towards the lowest energy state while sampling higher-energy states to ensure the global landscape is sampled.
The energies of the system after and before the MC trial are U new and U old , respectively. T is temperature and k B is Boltzmann's constant. The ensemble averages produce a picture of the system in equilibrium which is believed from extensive experience to be a unique thermodynamic equilibrium Algorithm 4.1 Iterative Scheme to solve the PNP-FS-IF system Given c 0,0 i , φ 0,0 i , for i = 1, · · · , N . for m = 0, · · · do for k = 0, · · · do Solve Nernst-Planck equation using EAFE for c m,k+1 Solve Poisson equation using standard FEM for φ m,k+ 1 2 , Update the electrostatic potential solution φ m,k+1 , φ m,k+1 = αφ m, end for Assign the solutions as initial data for the next time iteration, c m+1,0 end for independent of the initial conditions used to start the MC iteration. A Grand Canonical MC scheme [44] is used to maintain constant uV T , or constant chemical potential, volume, and temperature. It allows us to model a system of fixed volume and maintain a homogeneous bulk environment in the center of the system while adding more ions near the walls as the interfacial effects require.
Periodic boundary conditions are imposed on the directions parallel to the walls. This approach has been tested in previous work [65]. The chemical potentials needed to produce the desired bulk concentrations are determined in a separate Grand Canonical ensemble simulation [44].
The lengths of the cells were chosen to ensure that the bulk region amounted to at least half of the cell length. The cell width was adjusted to provide adequate sampling and computational efficiency. The number of MC moves for each simulation was 2.5 × 10 6 .
In the following subsection we present the computational results in uncharged wall and charged wall cases. 4.1. Uncharged wall cases. We consider an uncharged wall/boundary situation, that is, the permanent charge density is zero, i.e., ρ 0 = 0. Note that the electroneutrality constraint (24) has no effect in the bulk of the system (see the discussion in Remark 3). Thus, the function ψ should be constant everywhere in the domain. Moreover, the value of λ should be zero. Hence, the system (30), (31) is reduced to (18), (19). To maintain consistency and uniformity in our approach, one set of energy coefficients for the repulsive potential is chosen throughout the numerical experiments as ε i,j = 0.05 for i, j = 1, 2.
We compare the numerical results between PNP-FS-IF and MC computations. We consider the two ionic radii, a = 0.1, 0.2 (nm). The charge density results with monovalent and mono-divalent cases are presented in Figure 2, 3, respectively. PNP-FS-IF results show almost flat/constant behavior. At small radii, MC results show a desorption behavior near the wall due to surface tension, but as the ionic radius increases, packing effects increase the interface concentration above that of bulk. The difference may be caused by the diffusion effect induced by the Brownian motion in PNP-FS-IF but not in MC. The diffusion in PNP-FS-IF contributes to make the solution flat.
Remark 6. We need to investigate further the reason for the difference between PNP-FS-IF and MC even in uncharged wall situations. Moreover, it is not clear whether the solutions of the systems-PNP-FS-IF and MC-are a local or global minimum. In this case, a mathematical analysis and numerical computations are indispensable for pursuing the reasons for the differences. Discussion without analysis is unlikely to lead to definitive results. However, since the system includes a nonlocal and nonlinear term analysis is not easy.

4.2.
Charged wall cases. In this subsection, we consider the charged wall situations, i.e., the nonzero permanent charge cases. The charged wall is designed by the permanent charge, ρ 0 = 0.1(C/m 2 ) on both side walls. Thus, we compute the system equations, (30), (31) with the interfacial constraint. The function ψ is established by where d(x) is the distance function for the nearest wall, and η is the thickness of diffusive region which is related to the the thickness of double layer from the wall. The profile of the function ψ is presented in Figure 1. The value of λ is empirically chosen as λ = 0.001, which allows the electroneutrality constraint term to produce the charge inversion phenomena. We also consider monovalent-monovalent and mono-divalent salts, that is, z p = 1, z n = −1 and z p = 2, z n = −1. We compare the numerical results in cases of a = 0.1, 0.2 (nm). The monovalent cases are presented in Figure 4, and the mono-divalent cases in Figure 5. Both PNP-FS-IF and MC show a constant behavior in the bath region, and charge accumulation near charged walls. In these cases, the PNP-FS-IF shows steeper behavior near wall than the MC, Figure 4, Figure 5. PNP-FS-IF shows the charge inversion (layering) phenomena near the negatively charged wall in the mono-divalent cases, Figure 5.
As we discussed and showed in previous sections, the systems, PNP-FS-IF, MC are not identical. There is a reduction of dimensionality in numerical computations, and the difference in entropic contribution of ions. However, PNP-FS-IF is able to reproduce charge inversion phenomena similar to those in MC. 5. Discussion. We note that the comparisons between modified PNP and MC are made in reduced situations comparing three dimensional results to one dimension results. The new modified PNP is performed in one dimensional space, while, a cross sectional result is taken from three dimensional MC results. Thus, the systems are not identical and that may cause some difference between modified PNP and MC computations. Moreover, the PNP-FS-IF theory includes the Brownian motion of ions. That is, the total energy of the system in PNP-FS-IF theory includes the entropic energy, which produces the diffusion term in the resulting system of 'Euler Lagrange' equations. MC deals with the entropy in a quite different way. Finally, MC handles the electrostatics in a well precedented and justified way, but the charged sheet approximations of Torrie and Valleau [63] do not identically satisfy Poisson's equation for charged spheres with these boundary conditions. All of these facts may make a difference near the wall between PNP and MC in both charged and uncharged wall situations. 6. Conclusion. The charge inversion (layering) phenomena is a complicated and important behavior of ionic mixtures near a charged wall. It is difficult to recover the charge inversion in regular PNP theory. The energy variational approach provides a concrete process for the modeling of physical phenomena using a mathematical interpretation and energetic point-of-view. Using the energy variational approach, we recover the charge inversion phenomena in ionic mixtures using a new modified PNP equation (PNP-FS-IF) with a relaxed electroneutrality condition as an interfacial constraint. The new model recovers the charge inversion phenomena of ionic solution near charged walls in one dimensional space. Moreover, although PNP-FS-IF and MC computations are not identical, they show compatible results for the charge inversion (layering) phenomena. The PNP theory with the excluded volume effect needs further investigation in high dimensional spaces, but one dimensional models will always have important roles. High dimensional verifications will be reported once problems of computational efficiency and accuracy are solved.
It seems unlikely that the full range of layering phenomena important in chemistry, physics, and biology can be easily captured by any single theory. The range of concentrations, types and charges of ions, and the very large voltages used in electrochemical engineering imply that some phenomena will depend on both macroscopic and quantum chemical scales, along with the non-ideal properties of ionic solutions analyzed here. The variational method seems well suited to mix this diversity of physics. Relaxation of electroneutrality may prove a useful tool for understanding a much wider range of layering phenomena than we investigate here, including ion channels, batteries, and semiconductor devices.