Electric Double Layers with Surface Charge Regulation Using Density Functional Theory

Surprisingly, the local structure of electrolyte solutions in electric double layers is primarily determined by the solvent. This is initially unexpected as the solvent is usually a neutral species and not a subject to dominant Coulombic interactions. Part of the solvent dominance in determining the local structure is simply due to the much larger number of solvent molecules in a typical electrolyte solution.The dominant local packing of solvent then creates a space left for the charged species. Our classical density functional theory work demonstrates that the solvent structural effect strongly couples to the surface chemistry, which governs the charge and potential. In this article we address some outstanding questions relating double layer modeling. Firstly, we address the role of ion-ion correlations that go beyond mean field correlations. Secondly we consider the effects of a density dependent dielectric constant which is crucial in the description of a electrolyte-vapor interface.


Introduction
The occurrence of electrolyte solutions meeting solid (or liquid), materials is ubiquitous. The behavior of the resulting interfaces is key in the fields of surface and colloid science, electrochemistry and corrosion and soft materials and biomaterials. These interfaces are subject to both physical and chemical interactions with the fluid. Examples include chemical reactions between surface terminal groups and certain solution components, physical adsorption or desorption of charges at the interface [1,2]. All these interfacial phenomena may manifest themselves in the accumulation of surface charge at the interface. Consequently, the electrostatic surface potential will change. As a result, mobile charges in the vicinity will redistribute to lower the total free energy. Thus, surface chemistry effects propagate into the neighboring phase. The entire system, consisting of charged interface and locally distributed mobile charges and electrostatic potential is referred to as the electric double layer (EDL). The double layer name has historical roots, reflecting the initial notion that the charged interface with the electrolyte solution can be represented by a simple capacitor, where the surface charges and the solution counter-charges are placed on two well-defined planes in space [3]. This is certainly not an accurate representation of the physical situation, as entropic effects cause the charges to distribute themselves in space, away from a strict capacitor plane. An alternative description was suggested by both Gouy and Chapman [4][5][6] who realized that the EDL had diffuse aspects. They suggested a model based on Maxwell theory of electrodynamics, relating the potential Ψ(r) to the charge density distribution ρ e (r).
In parts of this paper we account for the presence of the solvent but we stress that this accounting is incomplete. It captures the important basic aspects of (1) excluded volume and (2) solvent-ion attractions (i.e., ion solvation) through the the Lennard-Jones (LJ) interactions. For instance, the effects of solvent orientation and solvent polarizability are not included in any detail, but rather mimicked through an effective LJ interaction. Similarly, the network structure of a solvent like water, say, including the quantum nature of hydrogen bonding are not captured at this time with current classical density functional theory (cDFT) techniques. Some authors have criticized models that do not include all known aspects of the solvent. We point out that there is always great value in the usage of simple models as their strength is to highlight and explore fundamental aspects and to provide insight into the phenomena of interest (e.g., the behavior of EDLs). In addition, there are recent approaches (see References [7][8][9]) that use experimentally measured ion solvation energies to improve the values used to approximate the LJ interactions, say. Such a strategy can address ion-specific effects that experiments have established.
The inclusion of more detail is then made when it becomes clear that the simple model approach fails. A nice illustration of such a situation can be found in the work of Wilson and Madden [10], who showed that to obtain the correct crystal structure for simple divalent salts like MX 2 it is necessary to include ion polarizability into the lowest free energy calculation.
Having recognized the strengths and weaknesses of density functional theory (DFT) it is worthwhile to stress that molecular simulation is the more convenient approach to explore the effects of solvation details including polarizability, orientation, hydrogen bonding and quantum effects.
The remainder of this paper is arranged as follows. In the next section we introduce the chemical boundary condition known as charge regulation and then in the next section we introduce the details of density functional theory that includes a neutral solvent component. We then discuss the role of screening ion-ion correlations that go beyond mean field. Finally, we address the role of the dielectric permittivity at an interface and introduce the Oleksy-Hansen method of a permittivity that depends on the local (coarse-grained) density. In the results section we illustrate the effects of charge regulation and the solvent inclusion by looking at the net charge distribution in a double layer and the corresponding electric potential distribution for three different wetting cases. We then illustrate the effects of including screening ion correlations for the primitive model and show how it can produce charge inversion. The latter is an impossibility in the Poisson-Boltzmann approximation. Finally, we illustrate the effects of ion correlations on fluid flow. The paper then concludes with a conclusions section.

Charge Regulation
Poisson's equation [11] is an example of a second-order partial differential equation (pde), When the right-hand side is equal to zero the equation is reduced to Laplace's equation. Poisson's equation describes the variation of the electrostatic potential, Ψ, in space, given a spatial distribution of charges, ρ e (r). In the study of electrolyte solutions the situation is slightly different, there we seek to find the equilibrium distribution of charges. This is a far more complex problem. Two familiar examples are the ionic distribution in a bulk electrolyte fluid and the distribution of ions in an electrical double layer. The first allows the determination of ionic activity coefficients whereas the second example leads to the prediction of forces acting between charged surfaces, as encountered in colloid stability studies.
To make progress with solving the Poisson equation one requires a simplification. This can be accomplished by postulating a relationship between the charge distribution and the electrostatic potential. Inspired by an ideal gas in an external field the approximation used stipulates that the density distribution is proportional to the Boltzmann factor of the potential. This approximation becomes exact in the low-density limit. Substituting this expression into the Poisson equation produces the so-called Poisson-Boltzmann (PB) equation, a nonlinear second order partial differential equation for the potential: where ρ 0 i is the number density of ionic species i and q i is the ionic charge including the sign (in units of e). The charge screening is characterized by the inverse Debye length κ is defined as, Debye and Huckel solved the linearized PB equation for an electrolyte solution whereas Gouy and Chapman solved the linearized PE equation for a charged surface. To solve the latter PB equation requires a boundary condition (bc). Natural and numerically convenient choices for bc are the potential at the surface (i.e., Dirichlet condition) or the spatial derivative of the potential at the surface (i.e., the surface charge, known as the Neumann condition).
In 1971, Parsegian and Ninham [12] proposed a radically different approach that is more true to the actual local chemistry. They drew attention to the fact that charged surfaces acquire their surface as a result of an ionizable surface group. That ionization is, ultimately, a chemical reaction where a neutral surface group splits of an ion leaving a charged surface group. A simple example would be a deprotonation of a surface group AH, viz.
This leaves a negatively charged surface. An equilibrium constant K (or equivalently a Gibbs free energy change of reaction, ∆G) relates the local concentrations of the reactants, here AH, A − and H + , where it is understood that both AH and A − are surface bound. In what follows, ∆G will typically be taken from the literature or set equal to the value in the bulk.
The relevant local concentration of the proton is that at the surface. If the surface is represented by a hard wall then this can be the contact density at the wall. On the other hand, if the wall is manifested as a smooth potential there is no contact density and a choice needs to be made. In previous work we have selected a density average over a narrow range, we shall employ that definition here. Alternatively, one could consider a profile based equivalent to the bulk y-function (i.e., y(r) = g(r) exp(φ(r)/kT)) [13], that is This is simply a smooth and continuous extension of the profile into the wall region. The contact density must then be measured at an effective hard wall position. The reader is referred to Reference [2] for more details regarding the use of the y-function in reactive systems and the use of y(z) in interfacial systems.
The Parsegian and Ninham approach replaces a physical bc (constant surface charge or constant surface potential) with a chemical bc, expressed in terms of a surface reaction's equilibrium constant (Ka or pKa). It is known as charge regulation (CR) as the surface charge is regulated by Ka (or pKa). This implies that both the surface charge and the surface potential are no longer input parameters but instead found from the solution of the electrostatic double layer problem. CR was originally proposed for a PB approach to double layers. More recently, Fleharty et al. [14] took the step to implement charge regulation in a DFT formulation of a double layer.

The Grand Thermodynamic Potential
EDLs involve multicomponent solutions, containing charged (ionic) and, if included, neutral (solvent) species. These species are subjected to the effects of the field exerted by the interface with the substrate. The field has electrostatic and van der Waals components, which lead to a variety of interactions with the fluid. In addition, the solution components (M) in the fluid interact with each other. All of these interactions are easily incorporated in a cDFT model. The electrolyte solution is described in terms of a grand thermodynamic potential functional, which for single flat EDL reads [14][15][16][17][18][19][20]] The first term on the right hand side of Equation (6) corresponds to the ideal contribution to the free energy, where λ = h 2 /(2πm i k B T) is the thermal de Broglie wavelength, h is Planck's constant, m i is the mass of species "i" and ρ i (z) is the local density of component "i" along the z coordinate normal to the wall. The radial coordinate R runs parallel to the EDL interface with the substrate. The excess free energy consists of hard-sphere and long-range parts. The hard-sphere contribution is based on the derivation of Rosenfeld [21,22] and reads Φ HS {n α (R, z)} is the hard-sphere reduced free energy and n α (R, z ) = n α (r) (r being the position vector) is the weighted local density.
Remarkably, the functional form of the reduced free energy Φ HS {n α (R, z)} is independent of the number components M. It is the weighted densities n α (r) exhibit such a dependence and are defined by and the weighting functions ω i α (r) are Here Θ and δ denote the Heaviside step function and Dirac delta function, respectively. The long-range contribution to the free energy functional is where Φ LR (R, |z − z |) is the long-range interactions contribution of the reduced free energy. Equation (12) indicates that the long range interactions are accounted for in a mean-field limit [23].
The term F SC [{ρ i (z)}] refers to the free energy due to ion-ion correlations and is discussed in detail below. The last term in Equation (6) is the Lagrangian constraint, which accounts for the external fields V ext i and fixed chemical potentials µ i for all species. The densities of all components ρ i (z) are found by minimizing the functional [15,24] which produces the so-called Euler-Lagrange equations. Typically, these are solved for the unknown component profiles, ρ i (z), by an iterative method. The charge density distribution in the EDL is derived by summing over the individual charged species, taking into account the charge numbers with their sign This approach includes the contributions of the interactions between the solutions species, which are present in the terms F ex [21,22], which accounts for the excluded volume effects and F ex long [{ρ i (z)}], which captures long-range interactions. In our model, these interactions are of the Lennard-Jones (LJ) [25] (15) and Coulombic type [14,[17][18][19][20], where d ij = (d i + d j )/2, d i is the diameter of component "i" and r ij (R, z) is the distance between species "i" and "j". All non-Coulombic interactions of a molecule (or ion) of type "i" with the interface are assumed to be of the "hard wall" type, that is, only the excluded volume effects are taken into account. The only exception is our analysis of the solvent-wall interactions presented in Section 3.3 below, where the LJ (9-3) potential [25] is used in the analysis. The electrostatic interaction of the ions with the interface is Other choices for the attractive interactions are also available. For example, Oleksy and Hansen [26][27][28] used a Yukawa potential to account for the long-range attraction between the solution species.

Higher Order Electrostatic Correlations
In the results presented so far, the electrostatic contribution to the energetic components in the cDFT formalism has been treated only at the mean field level; that is, no electrostatic correlations between the ions have been included except that they feel the average electrostatic potential. However, higher order electrostatic correlations beyond the mean field can be quite significant. This is especially true for multivalent ions but these correlations also play a role for monovalent ions. In this section, we review the origin and effect of these higher order electrostatic correlations.
To give an intuitive description of the origin of the electrostatic correlations, consider a homogeneous (bulk) electrolyte solution. There, the ions interact via the Coulomb potential, perhaps moderated by the dielectric constant of the solvent. Because the system is homogeneous, the mean electrostatic potential is identically zero; averaging over a sufficiently long period of time, there is charge neutrality everywhere. However, the nanoscale structure of how the ions arrange around each other (i.e., how the screen each other) depends strongly on the electrolyte concentration and the charges on the ions. Specifically, the length scale of screening decreases as either the concentration or the ion valences increase. This is reflected in the Debye length, the most commonly used estimate of the screening length. If ions are generally closer to each other when they are at high concentration or when one species is multivalent and this is not reflected in the mean electrostatic potential, then there must be an energy term in the cDFT formalism (Equation (6)) that accounts for these correlations between the ions. We call this the screening term, denoted with SC.
Virtually all modern cDFT formulations of the screening energy term are based on the Mean Spherical Approximation (MSA). The MSA is a theory of homogeneous electrolytes that extends the classical Debye-Hückel theory to include the size of the ions [29][30][31]. Other theories like the Hypernetted Chain approach are more accurate [31] but the MSA has the advantage of explicit analytic formulas, which makes it particularly convenient to generalize it to inhomogeneous electrolytes in cDFT.
Here, we do not focus on different cDFT formulations that have been created to approximate this term but rather the effect the screening term has on EDL structure. A number of different expressions have been derived for the screening term over the last 30 years. These provide different levels of trade-off between accuracy (EDL structure compared to Monte Carlo simulations) and computational speed. For the results shown here, we use the Reference Fluid Density functional of Gillespie et al. [32,33], which remains the most accurate one to date. Other cDFT formulations include the bulk reference Taylor expansion [34,35] and the functionalized MSA [36]. The relative accuracies of these three approaches was recently catalogued [37].

Molecular Interactions, Solvent Polarity Effects and Dielectric Permittivity
Electric double layers typically form at interface of substrates with electrolyte solutions. Hence, the Coulombic interactions that involve the ionic species are of primary importance. However, the existence of a stable liquid phase requires the presence of attractive forces forces between the solvent molecules. The solvent molecules also interact with ions, which accounts for the ionic solvation. Depending on the solvent polarity, its molecules may be involved in isotropic attractive (e.g., Lennard-Jones, Yukawa, etc.) or dipole interactions. In addition, all species have finite size and that leads to a short-ranged repulsive force, which plays a major role for the liquid structure at the molecular scale.
A simple description of electrolytes and EDLs can be accomplished by ignoring all interactions except for the Coulombic (and in some cases the short-range repulsion) between the ionic species. Such models are colloquially referred to as "primitive" [38]. The primitive models [36,[39][40][41][42][43][44] consider the solvent to be a structureless continuum, characterized by a uniform relative dielectric permittivity ε r . The model reflects the ionic interaction contributions but fails to capture the structural effects caused by the presence of the neutral solvent. At the opposite extreme are the "civilized" models [45][46][47][48][49][50][51][52][53] that account for the presence of all solutions components and all possible interactions between them-Coulombic, van der Waals, dipole-dipole, ion-dipole and so forth. The inclusion of the dipole effects, that are due to the solvent polarity, requires an additional integration to properly average all dipole orientations. Such a step needs significantly greater computational resources and time. A reasonable compromise is offered by the "semi-primitive" models, which explicitly account for the presence of solvent molecules but ignore the effects due to the their polarity. This approximation means that the solvent molecules exhibit short-range repulsive forces and isotropic attractions (necessary to ensure the existence of a stable liquid phase) but are not involved in any orientation-dependent dipole-dipole or ion-dipole interactions. The semi-primitive models proved themselves helpful in demonstrating the effect of the liquid structure on the properties of charged electrolyte interfaces [14,[17][18][19][20][54][55][56][57][58][59]. Since the semi-primitive approach does not account for any dipole effects, the dielectric permittivity has to be added ad hoc in order to properly scale all Coulombic terms. A very insightful analysis of the problem was recently offered by Oleksy and Hansen [26,27], who argued that this approximation is quite reasonable and can be further improved by introducing a dielectric permittivity that depends on the local weighted solvent densityρ 0 (z) defined bỹ The local weighted solvent density density can then be used to calculate the local dielectric permittivity ε r (z). An example is the Clausius-Mossotti equation where m is the molecular dipole moment. The Clausius-Mossotti equation can be derived using a mean-field cDFT model of dipolar hard sphere fluid [28]. Unfortunately it fails for high dipole moments and cannot be used for liquids such as water [27,28]. For such cases, Oleksy and Hansen [26,27] proposed the empirical expression The function ρ 0m (T) = [ρ g 0 (T) + ρ l 0 (T)]/2 is the mid-point fluid density between coexisting gas and liquid phases [with densities ρ g 0 (T) and ρ l 0 (T), respectively], whileρ 0 (z) is given by Equation (19). The empirical function f (T) = 88 − 0.37T, was designed to conform to experimental data for water permittivity.
The Oleksy-Hansen approach is particularly useful when the semi-primitive cDFT model is applied to wetting electrolyte liquid films in coexistence with a gas phase. In such cases, the dielectric permittivities in the liquid and gas phases may vary by more than an order of magnitude and that is adequately captured by Equation (21) above. This model was tested against the civilized cDFT model proposed by Biben et al. [52], which fully accounted for the dipole effects at the molecular level and showed a very good agreement. Single phase liquids are less challenging since the dielectric permittivity does not experience such abrupt changes. Still, the local dielectric permittivity in a liquid solution is perturbed by the presence of charge (i.e., an ion or charged group.) This perturbation decays with distance and after a few molecular diameters the permittivity resumes its bulk value [51,52]. Hence, the detailed solvent polarity effect on the local dielectric constant may be neglected for low to moderately concentrated electrolyte solutions. For example each two ions in a 0.01 M solution are separated, on the average, by more than 30 solvent molecules, which is more than sufficient distance for the permittivity to relax to its bulk value. At higher electrolyte concentrations and/or in the presence of multivalent ions the local variations in the local dielectric permittivity should be taken into account either using the Olesky-Hansen approach (see Equation (21)) or by developing a full-scale civilized model including all possible interactions (van der Waals, Coulombic, dipolar, etc.). Another feature of concentrated electrolyte solutions is that the ionic screening correlations become important and need to be properly taken into consideration [36,44].

Solvent Effects
As an illustration of the importance of including the solvent in cDFT we consider an EDL in a situation where the solvent is attracted to the wall with three different wettability conditions. Thus, we mimick a solvophilic wall, a solvophobic wall and a partially wetting wall that we will refer to as a neutral wall. The wetting variations are accomplished by varying the strength of the LJ interactions (i.e., through ε s−w -see Equation (17)) between the wall and the solvent component. We vary ε s−w /k B T from 0 to 1 to 2. Note, that for the ion species the value is held constant at 1 throughout. Although we can be sure that the case ε s−w /k B T = 0 produces complete drying at liquid-vapor coexistence, the other two values are merely inspired guesses and we are not implying that (at liquid-vapor coexistence) they would produce complete wetting for ε s−w /k B T = 2 or a contact angle of 90 degrees for the neutral wall. The results for the spatial net charge distributions are shown in Figure 1. These are the results for a symmetric monovalent electrolyte at fairly low molarity (0.01 M). The bulk pH is set at 4 and the charge regulation parameters are set at pK + = −2 and pK − = 6. The solvophilic (red) curve, ε s−w /k B T = 2, shows the most pronounced build-up of positive charge variation in the EDL. Note that all three wetting states display a pronounced layering oscillation with a period of 1 solvent diameter. As mentioned this is due to the packing of the solvent molecules. The increased number of solvent molecules at the wall prevents the PDIs form approaching and neutralizing the negative surface charge. Hence, the more negative interface attracts the positive counterions in the EDL. This layering modulation is carried over to the positive and negative ions and hence is reflected in the profiles of the net charge (ρ e (z)).
The corresponding profiles for the electrostatic potential is shown in Figure 2. As expected the largest value for the surface potential (i.e., the potential value at z = 0) occurs for the solvophilic case, ε s−w /k B T = 2. It is interesting to see that the electrostatic potential is monotonic throughout the EDL. That is, there is there are no hints from the layering seen in Figure 1. This is consistent with Poisson's equation (Equation (1)), which shows that it is the second derivative of the potential that corresponds to the (net) charge distribution. Since the overall sign of the fluid charge density does not change, the potential curvature (or d 2 Ψ/dz 2 ) will not change as well. The magnitude of the curvature, however, will change according to the the fluid charge density curves shown in Figure 1 but this is hard to notice by visually observing the potential curves in Figure 2. However, differentiating the potential curves twice will recover the oscillating results for the fluid charge density ρ e (z) shown in Figure 1. The results shown here are a demonstration of the role of including an explicit solvent in cDFT. It is clear that the solvent's presence is necessary to capture the effects of different wettabilities of the solid phase. Moreover, to capture thin wetting films it is essential to include the formulation of a density dependent dielectric permittivity described in the previous section.

Effect of Screening Correlations on EDL Structure within the Primitive Model
The vast majority of the studies studying the effect of the screening term were performed with the primitive model of ions, where the ions are charged hard spheres and the solvent is a background dielectric material; there are no explicit water particles, as there are for the other results shown in this article. Therefore, it is not known exactly how the inclusion of higher order electrostatic correlations will affect these results but a number of generalities will carry over. The focus of this section will be discussing these general properties but within the framework of the primitive model of electrolytes. While many of these properties have been generally known for a while, how they all come together to define EDL structure was only recently described in detail by Voukadinova and Gillespie [44].
In the primitive model, ion concentration profiles as a function of the distance x from a smooth, hard, uniformly charged (as opposed to the charge regulated approach described above) surface are given by where i is the ion species. Here, we assume that the surface charge is negative and constant (i.e., not regulated). The three terms in the Boltzmann factor are different components of the interaction physics. The first term (HS) is the hard-sphere contribution, and ∆µ HS (x) is the energy it takes to insert an uncharged ion at distance x, relative to the bulk value far from the wall. Similarly, ∆Ψ (x) is the mean electrostatic potential relative to bulk. The third term is the screening term. All these terms are the functional derivative of a corresponding cDFT free energy term, see Equations (6) and (13).
To see the effect of the screening term most clearly, we rewrite this to focus on the electrostatic potential profile. Specifically, we write this in terms of the counterion (cation) species but for convenience drop the species notation: In this equation, the first two terms on the right-hand side are negative. The first term is negative because, when the surface charge on the wall is large enough, the cation concentration will be greater than the bulk value far from the wall. This also causes the second (HS) term to be negative because the more dense an area is, the more difficult is to insert an uncharged particle. The SC term, on the other hand, is almost always negative [44] because, just like in a homogeneous system, having higher ion concentrations makes the screening energy more negative. This makes the third term positive.
In the classical Poisson-Boltzmann theory where the HS and SC terms are ignored, the electrostatic potential given just by the first term in Equation (23). Therefore, the balance between the HS and SC terms determines how the electrostatic potential is different from the classical case. In their recent work, Voukadinova and Gillespie [44] found that at low to moderate service charges the SC term in Equation (23) is more positive than the HS term is negative; the balance is only flipped at extremely high surface charges when the counterion concentration is so large that it becomes extremely difficult to find space for more ions. Therefore, the effect of the higher order electrostatic correlations is to make the electrostatic potential less negative than the classical Poisson-Boltzmann case.
In fact, the screening term can become so large as to make the electrostatic potential positive. This is known as charge inversion. Charge inversion occurs mostly for multivalent ions and this is because the screening term scales as the square of the counterion valence [36]. Thus, the screening term becomes significantly larger the higher the valence.
The effect of the screening term in making the electrostatic potential less negative than the Poisson-Boltzmann case and leading to charge inversion is shown in Figure 3a. There, the onset of charge inversion is shown by increasing the bulk concentration of an electrolyte with divalent cations. At low concentrations (blue lines), the electrostatic potential with higher order electrostatic correlations (solid lines) is less negative compared to Poisson-Boltzmann without these correlations (dashed lines). This trend continues at high concentrations (black lines) but the potential with screening correlations becomes positive (solid line). Because the Poisson-Boltzmann theory does not include this negative contribution it can never have a change of sign in the electric potential.
This positive electrostatic potential has two consequences. Right now we focus on the consequence for the EDL structure; below we consider the effect on fluid flow. Specifically, the positive potential draws in anions (co-ions) in a second layer of ions behind the initial high concentration layer of cations (counterions). This is shown in Figure 3b. There, the anion concentration (solid red line) increases above the bulk concentration (thin red line); for the Poisson-Boltzmann case (dashed red line) this does not happen. This negative-positive-negative sandwich of surface charge, counterions and co-ions is a hallmark of charge inversion and it is due to the change of sign of the electrostatic potential caused by the higher order electrostatic screening correlations. This is the most extreme consequence of these correlations but in general they make the electrostatic potential less negative and therefore tend to decrease the cation concentration and increase the anion concentration, relative to the classical uncorrelated Poisson-Boltzmann theory.

Effect on Fluid Flow
Beyond this effect of the electrostatic screening correlations on EDL structure, they and the concomitant charge inversion can also have a significant effect on fluid flow. This was explored in detail in References [60,61] for nanofluidic channels where electrolytes flow parallel to two charge surfaces that are of the order of 100 nm apart [62]. The geometry is illustrated in Figure 4  In such a system, advection moves water through the channel based on the local mean electrostatic potential. (This movement of fluid is in addition to the movement of ions.) Specifically, the advective velocity profile at each longitudinal location y across the channel that drives fluid down the axial direction x of the channel is proportional to Ψ (y) − Ψ (H S ) where H S is the slip plane location where the velocity is zero [63]. If one considers the velocity at the center of the channel where, for convenience, we take the electrostatic potential to be zero, then the direction of the velocity is given by the sign of −Ψ (H S ). Consequently, if there is charge inversion and the slip plane is in the region a positive potential, then the fluid will flow in the opposite direction of when there is no charge inversion and also in the opposite direction of the ions.
Each of these has a consequence that can be measured in experiments. In the classical case without charge inversion, the counterions in the EDL and the water move in the same direction when a hydrostatic pressure gradient or electrostatic potential gradient is applied down the axial x direction of the channel. However, when charge inversion moves the fluid in the opposite direction of the ions, this reduces the total ion current (as some ions are in each volume of fluid that is moving) and can even change the sign of the current. These were the first measurable consequences of charge inversion in nanofluidics [64,65] and can be reproduced by cDFT [66,67].
In more recent work [60], the two wells at the end of the channel contained different electrolytes, one that exhibited charge inversion and one that did not. When an electrostatic potential was applied between the two wells that moves the two fluids towards each other (as the charge inverting fluid moves in the opposite direction of the other fluid), a stable front between the two fluids is established. Moreover, at this junction, low-concentration analyte ions will accumulate so that they can be pre-concentrated for analysis later. While the explanation of the physics of both the stable front and the ion accumulation is beyond the scope of this paper (see References [60,61]), this practical application is a macroscopic manifestation of the charge inversion produced by the higher order electrostatic screening correlations.

Conclusions
EDLs are inherently complex interfacial structures involving a multicomponent fluid with chemical ionization reactions present at the surface. To gain insight into EDLs it is often helpful to also consider basic models, as we have set out to do in this paper. One important aspect concerns the electrostatic screening effects acting between charges. These have been successfully studied within the primitive model, where the solvent is structureless background fluid. However, the conclusions reached by using the primitive model are transferable. That is, the main effect that the electrostatic screening correlations make the electrostatic potential less negative will still apply when explicit water particles are included. Moreover, they will remain an important contributor to charge inversion. However, what the relative balance of the screening term and the other terms will be when water is included remains to be determined. Nanofluidic flow was presented an application of the relevance of charge inversion.
EDLs with charge regulation boundary conditions and an electrolyte that contains an explicit solvent component display different rich interfacial behavior. That behavior includes a detailed structure dominated by the neutral solvent molecules, that strongly influences where the charged ions will be positioned. The solution structure near the charged interface affects the surface chemistry and hence there is a strong coupling between the solvent-induced effects and the charge regulation (for more details see References [14,18]). As an example of solvent effects we have presented the results for three types of wetting situations: solvophilic, neutral and solvophobic. These wetting types were created by enacting different interaction strengths between the solvent and the wall (while keeping all other interaction fixed). We find that although only the interactions between the wall and the neutral solvent were varied, the charge and electric potential were directly affected.
In this paper we have exclusively focused on simple planar surfaces. However, we stress that the methods discussed here are not limited to these simple cases. In fact, they can and have been applied to biological problems such as as ion channels. An example is the work of Gillespie et al. on the ryanodine receptor ion channel, see for instance, Reference [68].

Acknowledgments:
The authors are grateful for helpful discussions with Raviteja Vangara, and would like to thank Reviewer 2 for drawing our attention to References [7][8][9][10].

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript:

EDL
Electrostatic double layer cDFT Classical Density functional theory HS Hard sphere LJ Lennard-Jones MSA Mean spherical approximation HNC Hypernetted chain SC Screening correlations