The Dramatic Effect of Water Structure on Hydration Forces and the Electrical Double Layer

: Forces between hydrophilic surfaces mediated by water are important in various systems from lipid membranes and solid surfaces to colloids and macromolecules, first discovered as a significant addition to DLVO forces at the nanoscale. These “hydration forces” have been studied in great detail experimentally using osmotic stress measurements, surface force apparatus, and AFM, and they have also been the subject of multiple theories and simulations. One spectacular feature observed in experimental and simulation studies was the nonmonotonic, oscillatory decay in the forces between atomically smooth surfaces. Forces between “rougher” surfaces exhibit only quasi-exponential, monotonic decay. Here we revisit this hydration force problem by exploring the consequences of an extended phenomenological Landau − Ginzburg approach that describes nonlocal correlations in water, linking them with the key features of the wave-number k -dependent nonlocal dielectric function of water. With corresponding boundary conditions, this theory predicts the observed oscillatory decay in hydration force between ideally flat surfaces, the oscillatory mode disappearing with just a tiny roughness of the surfaces (of mean height ca. of the size of a water molecule). This study also brings an important side message. Explanation of these observations appears only possible under an assumption of two modes of polarization in water, consistent with the behavior of the response function, i.e., Lorentzian at small k and resonance-like at higher k . This resolves the “force oscillation − non-oscillation” paradigm, which is a strong, although indirect indication of the existence of these two modes. We also consider other important subjects, such as how the distribution of ions near a charged surface reacts to the propensity for overscreening oscillations due to polarized water. This is important not only for the interactions between charged surfaces but also for the fundamental understanding of the structure of the electrical double layer at electrochemical interfaces. We show that even in dilute electrolytes, the distribution of ions in the vicinity of the polarized interface follows, although not literally, preferential positions corresponding to the potential wells caused by “resonance” water layering. For a sharp interface, the theory predicts that the decaying spatial oscillation profiles extend over a 1 to 2 nm distance from the interface. With the smearing of the interface and the corresponding suppression of the resonance water layering, oscillations in the spatial distribution of ions subside, resulting in a familiar Gouy − Chapman − Stern picture. At longer distances from the interface, whether smeared or not, the ion distribution profiles become Gouy − Chapman-like. The effect of the boundary conditions on water polarization at the interface goes beyond a trivial shift of the potential of zero charge. We show that they can dramatically affect the ion distribution near the charged surface. Last, but not least, we study how the interfacial water layering influences the double layer capacitance and show the effect of the boundary conditions on the slopes of Parsons − Zobel plots, resolving some recently discussed puzzles.


I. INTRODUCTION
Hydration forces have had a long and dramatic history since their first experimental discovery−for a review, see the classical monographic account in ref 1, with some updates in ref 2. These forces have been experimentally studied between soft and solid surfaces 1 and even between macromolecules such as polysaccharides, collagen, and DNA. 3,4 Hydration forces, particularly between neutral lipid membranes, have received a lot of attention from biophysicists. 5 The main method in characterizing hydration forces between lipid lamellae was the measurement of osmotic stress isotherms. 3 Between solid surfaces, these force measurements were initially based on Surface Force Apparatus (SFA), 6 and later, Atomic Force Microscopes (AFM) 7 were also used to measure forces rather than just for drawing surface landscapes (see Section IV.B for a discussion of the latest applications of AFM techniques). The general understanding of the origin of hydration forces was that hydrophilic surfaces perturb the structure of water or affect its polarization at the boundaries and that these perturbations propagate from one surface to another through the correlations of the perturbed quantities in water. Thus, beyond the typical correlation lengths, such forces should vanish. Since many of the studied surfaces are charged, hydration forces are usually perceived as an additional force emerging at short distances on top of the forces described by DLVO theory−the force determined by the repulsion of electrical double layers (the socalled electrostatic force) and van der Waals attractions. 8 The mathematical description of such a liquid-mediated force was pioneered by Marcělja and Radic. 9 They postulated a simple gradient expansion of a free energy functional of the liquid (water), quadratic in the "order parameter" (often interpreted as the orientation of the molecules), where the coefficients determine the correlation length. Under the assumption of opposite orientations of the order parameter at the boundaries, such a description gives a repulsive force, exponentially decaying with the distance between surfaces with decay length equal to the correlation length. Cevc and Marsh 10 addressed attention onto the possibility of introducing softer boundary conditions into the Marcělja-Radićtheory; much later Cevc and Kirschner 11 explored the role of smearing of the boundary conditions on the predicted and measured forces.
Belaya, Feigelman, and Levadny (BLF) 12 related the description of hydration forces with the nonlocal dielectric response of water, 13,14 using a form of the boundary condition on the polarization of water on the surfaces and the Lorentzian form of the nonlocal, wavevector dependent dielectric function, earlier proposed by Kornyshev,13,14 the form being equivalent to the gradient expansion functional of Marcělja and Radic. A similar approach was also developed by Cevc. 15 Kornyshev 16 proposed a more general scheme relating the hydration force with the nonlocal dielectric function: losing its translational invariance at the interface, the problem is reduced to the solution of integral equations for the inverse dielectric function. He demonstrated this method on the same model of the bulk dielectric function (as used by BLF) and obtained qualitatively similar results, but more generic due to a more general form of the boundary conditions and the allowance for distortions of the nonlocal dielectric response near the interfaces.
Later the interest was shifted to account for the possible inhomogeneity of the surface boundary conditions, particularly important for the hydration forces between lipid lamellae with zwitterionic polar heads. Kornyshev and Leikin 17,18 developed a theory showing that the decay length of the propagation of the force is determined not only by the spectrum of correlations in water but also by their interplay with the correlations of charge density of zwitterionic groups along the membrane surfaces. This explained for the first time how the decay length of the hydration force could vary when acting between surfaces of differing lateral structures. Follow-up work accounted for the effects associated with the rearrangement of the soft boundary conditions in response to the hydration interactions. 19−21 In the 1990s, it became clear that for water the Lorentzian approximation of the static nonlocal wavenumber dependent dielectric function was not sufficient (see ref 22 and the references cited there). Cherepanov 23 linked the hydration force with a more sophisticated form of the nonlocal dielectric function of water, having borrowed and fitted its form from the results of Bopp, Kornyshev, and Sutmann (BKS). 24 He showed that including this overscreening resonance in the domain of wavenumbers about 2π/d where d is the diameter of the water molecule gives rise to the decaying oscillation form of the hydration force.
In the "primitive" gradient expansions representing the Landau free-energy functional of bulk water, the hydration forces were always exponentially decaying with the distance. However, experimental crossed-cylinder SFA measurements of the forces between ideally smooth mica surfaces have shown sharp decaying oscillations, with a period on the order of the diameter of a water molecule. 25, 26 These oscillating patterns clearly have to be described by going beyond the simple gradient expansion of the free energy, or by using resonance-containing forms of the nonlocal dielectric function. 23 However, for atomically rough surfaces, these oscillations disappear. 27,28 Including the resonance alone does not explain why for some interacting surfaces the forces in the experiments decay exponentially, whereas for many others we see these oscillations. Wisely, it was suggested that smearing of the interfaces due to surface roughness leads to oscillation dysphasia, eliminating these oscillations. The question remained open, how large must the smearing be for this, and what would remain of the force after such smearing?
Also, the role of electrolyte ions needed better understanding. It is well-known that their presence gives the renowned electrostatic contribution to DLVO theory. The repulsion of electrical double layers 6,8 extends over the Debye screening length, which can be much longer for diluted electrolytes than the extent of action for hydration forces (for a 1 mM solution, the Debye length in water is ∼9.6 nm). Thus, the hydration force manifests itself only at short, nanometer distances on top of that repulsion. Often the experimental data for hydration forces were presented as a result of subtracting the theoretical DLVO force from the measured "total" force, 8 decoupling the two contributions. Intuitively, such a procedure is likely to be valid at low electrolyte concentrations when the Debye length is much larger than characteristic lengths of the hydration force. However, is it right to treat them as decoupled, or can the electrostatic interactions of the overlapping electrical double layers somehow affect the features of the hydration interaction?
It is natural to expect such an influence if the concentration of ions, or the charge of the surfaces, is high. Early theoretical works that predicted and described the observed monotonic decay of the hydration force in the linear response approximation have found that the presence of electrolyte can modify the decay length of hydration forces, 29 as well as the order parameter correlations affecting the Debye length. 30 Even more complex implementations of that interplay have been envisioned under conflicting conditions on surface charges and the boundary values of the order parameter. 31 However, none of those predicted effects, weak and strong at low and high electrolyte concentrations, respectively, have been confirmed experimentally. 1 Thinking beyond simple exponential decay, if the force is oscillating in space due to overscreening, would it not be reasonable to expect cations and anions to get trapped in corresponding extrema of the potential, which is also oscillating in space? How would the accommodation of ions there affect the resulting potential and the force between the surfaces? Would it increase or decrease the latter, amplifying or suppressing the amplitude of the force oscillations? Reference 23 has reported some results on this point, and we will discuss these toward the end of this paper. However, the primary question remains�is it possible to decouple the hydration interaction from the interactions of the overlapping double layers and, if "yes", under which circumstances? Indeed, as in many reported measurements the extracted profiles of hydration forces do not seem to depend on the concentration and type of electrolyte ions, unless there is a substantial difference in their propensity for specific adsorption onto the surfaces. 1 This paper is devoted to answering these and several other questions. Altogether, we summarize them as follows. Some of these questions have been recently addressed by De Souza et al. 33 Blossey and Podgornik have also recently developed a general field-theoretical approach to this problem, exploring this ion charge density and solvent polarization coupling, 34 and have studied various dipolar models of the liquid. Our approach below will take similar lines, but we will focus specifically on a phenomenological free energy functional that reproduces the simulated dielectric response of water. To forewarn the readers: in the phenomenological theory presented below, we will not go into the full complexity of the problem, which may require taking into account the lateral inhomogeneity of the boundary conditions, 1,17,18 fluctuations of the boundaries (potentially important for soft surfaces; for a review of a series of works on this subject, see ref 18), possible force-induced restructuring of the boundaries, 19,20 and the nonlinear response of water polarization. 35 We will nevertheless try to build a platform for comprehensive answers to the questions posed above.
Wherever possible, the results of the theory will also be compared with existing experimental data. This comparison will give us additional information for the complete picture of hydration interactions through layers of water, as well as of the structure of the electrical double layer at electrochemical interfaces. A reader less interested in the mathematical formalism with which we start below may proceed directly to the results of the theory in Sections III and IV.

II.A. Free Energy Functional for Polarization Fluctuations.
Below we formulate the principles of a phenomenological description of water confined between two polarized surfaces, which will be extended in Section IV to the description of the double layers forming at charged surfaces separated by electrolytes. The reader should not perceive this formalism as a molecular theory, but rather as a continuum field-theory platform which will allow us to investigate several expected but, most importantly, less expected consequences of the molecular correlations in water, as well as the dramatic role of the "conditions" on the surfaces.
We can construct a phenomenological energy functional describing the dielectric properties of a nanoconfined slab of water from the electrostatic energy depending on the displacement field D, and a correlation term written in terms of the polarization density field P tot : The vacuum electrostatic energy term in the free energy functional can be derived from the energy density of an electrostatic field, h E e 1 2 0 2 = | | , and the definition of the displacement field, D 0 = ε 0 E + P tot : where ε 0 is the vacuum dielectric permittivity. Approaching the correlation term in the free energy functional, we expand on the Landau−Ginzburg approach (LG) taken by Monet et al., 36 including higher gradient terms in the polarization density in the expansion. In the spirit of ref 37, we assume that the local polarization P tot is a sum of three contributions, P tot = P 1 + P 2 + P * , where P 1 and P 2 describe contributions from the reorientation of water molecules, and P * describes contributions from the polarizability of the molecules themselves. By introducing a three-component model, we can account for the different possible modes and spatial scales of polarization fluctuations. The intramolecular modes (P * ) include contributions such as fluctuations of electron density and intramolecular vibrations. However, we assume that these modes give purely local polarizability, independent of the wavevector, k, and so we do not include them in our model free energy functional. The other orientational modes could correspond to reorientation of molecules that results in rearrangements of the hydrogenbonding network, or oscillations about an average position preserving the hydrogen bond network. We avoid specifically assigning these modes here, and instead consider a model described by two different orientational modes P 1 and P 2 , one of which (P 1 ) we expand to the second derivative term in the LG expansion, and the other (P 2 ) to just the first derivative as in the Marcělja-Radic̀model, also known as the "Lorentzian approximation", where the set of phenomenological model parameters is given by {K ij }. As structural transitions at the surface are beyond the scope of this work, we keep the expansion quadratic in the polarization density. However, inclusion of the higher gradient terms in the LG expansion, as well as this two-polarization model, will allow us to build up a kernel that mimics the simulated nonlocal dielectric properties of bulk water. This model of two orientational polarization modes, correlated differently in space, one of higher frequency with decaying oscillations and the other one of lower frequency that decays exponentially, is an approximation. It matches some simulation data 22 as well as several consequences of the nonlocal dielectric response of water. 22,38 This concept has been discussed in detail in ref 22 in the context of the interplay between temporal and spatial correlations of polarization fluctuations. It is tempting to relate this splitting of the polarization density with the well-known attempts to fit the Debye spectrum of the frequency dependent dielectric permittivity of water with two different relaxation times, 39−43 one fast and one slow, often associated with librations of individual molecules and reorientations of their larger clusters, respectively. But this does not mean that water should be literally considered as a coupled system of two liquids, each with their own relaxation processes, although new arguments to validate such a picture have been raised 44,45 in the context of recently discussed "low" and "high" density clusters in water. But all these should better be interpreted in the language of collective modes. 22,46,47 II.B. Linear Response Functions. It should be noted that in this work, we do not consider dynamics of the response functions, that is, the derived response functions are static: i.e., ω = 0. For small electric fields, the polarization density P of the medium is related linearly to the displacement field D. As the response may be nonlocal, we can generalize the expression to P α (r) = ∫ V′ dr′χ αβ (r, r′)D β (r′), where α, β = x, y, z, and χ αβ is the dielectric response tensor. In a homogeneous medium, i.e., χ αβ (r, r′) = χ αβ (r − r′), we obtain an algebraic expression for the Fourier transforms, For isotropic media, we can decompose the response tensor χ αβ (k) into longitudinal (χ || ) and transverse (χ ⊥ ) parts with respect to the wavevector k; here we focus on the longitudinal component, for more details see ref 22. Hence, to obtain the K ij parameters of the expansion presented above, we can derive an expression for the static longitudinal dielectric response function from the free energy functional and fit it to simulated results. Variational minimization of the Hamiltonian leads to the following pair of coupled equations in Fourier space: Writing these equations in matrix form, and using eq 4, we obtain for the response function:  (9) We can hence also obtain the static nonlocal dielectric function ε(k), via the relation that χ || (k) = 1 − 1/ε(k): 13 It is important to note that the form of this expression for ε(k), and hence χ || (k), has a potential defect; it is not correct for any arbitrary set of parameters {K ij }. Depending on the parameter set chosen, ε(k) can violate the Dolgov−Kirzhnitz−Maksimov (DKM) rule; i.e., it can enter the forbidden band of 0 ≤ ε(k) < 1, which arises from stability arguments of the system. 48 This is one of the shortcomings of introducing a "minimal number" of terms in our phenomenological free energy functional. Hence, we must make sure that for the chosen parameter sets, ε(k) does not enter this region. However, we do not encounter this problem if we fit {K ij } to reproduce simulations of ε(k), the results of which do not violate the DKM rule. To parametrize this model, we fit the results of MD simulations of the nonlocal dielectric response of water with the broadly used SPC/E model for water (see Figure 1). Note that other rigid models such as Tip-np give similar results. The numerical values of the parameter set (ϵ w , k min , χ min , k max , δ k ) which determine {K ij } are found by focusing on three key properties of the nonlocal dielectric response of water, (i) its bulk permittivity, (ii) the position and value of χ(k) at low k, and (iii) the pseudoresonant peak of χ(k), indicative of overscreening. For part iii, we choose to fit the position and the peak in χ(k), which has been shown to capture quantitatively the range of polarization correlations. 36 This gives us five conditions on the parameters: Details for the simulation runs and data exploitation are provided in Appendix A.
We have neglected to include the "shoulder" (satellite peak) that is seen on the right side of the main peak, at k ≈ 5 Å −1 . (see also refs 22 and 24). This shoulder is a signature of another stationary eigenmode corresponding to the contribution of charge oscillations with a shorter period than the peak at k ≈ 3 Å −1 . This is likely to be associated with intramolecular charge structure. We can in theory incorporate this shoulder into the LG expansion at the cost of more terms, and hence more parameters to fit. In an unpublished work, we have considered a free energy functional expanded up to the 10th polarization derivative. Such an expansion can capture the values and positions of the main peak and the shoulder. However, the result of including the shoulder only affects the correlation lengths present in the system in a minor way, as it convolutes the oscillating envelope with a signal of shorter period and small amplitudes. For the purposes of this paper, however, we do not want to complicate the analysis, and we focus instead on the effect of the behavior at small wavenumbers and of the main peak. Additionally, since this shoulder occurs at higher wavenumbers, its effect would disappear first in the presence of even smaller roughness.
If we interpret the two orientational contributions (P 1 and P 2 ) to the polarization density as different "frequency" polarization modes in the liquid, then we can define a certain ε * , the dielectric constant at the frequency separating these modes. For example, in water, high frequency modes are attributed to librations and intramolecular vibrations, while the lower frequency modes correspond to molecular reorientations. The individual weighting factors of each contribution are encoded in the set of phenomenological parameters {K ij } which we have fitted to the simulation data. Hence, with this parameter set, and by finding the poles of the response function in (9), we can find the different characteristic lengths in the system. Here we find one pure decay length (λ d,2 = 3.65 Å) related to P 2 , and one decay length (λ d,1 = 1.71 Å) coupled with an oscillation ( λ o = 2.10 Å) related to P 1 .
It is however not clear in eq 9 what the relative contributions of the two polarization modes are. Hence, we need to express the response function in the form: where the weights a 1 and a 2 are related to the Pekar-like parameters: 49 a a where ε w = 71 is the macroscopic dielectric constant for SPC/E water. We can decompose eq 9 into the following form: where Q = k 2 , the coefficients c 10 , c 11 , and c 2 are dependent on the set of model parameters {K ij }, and Q 10 , Q 11 , and Q 2 are related to the roots of eq 9. The second term in the brackets of eq 13 can be attributed entirely to P 2 , and hence Q 2 = k d,2 2 = 1/λ d,2 2 . Therefore, the first term is the contribution entirely from P 1 . When Q → 0, χ ∥ (Q = 0) = 1 − 1/ε w = 1/a 1 + 1/a 2 . We can then see that Calculating these parameters for our fitted {K ij }, we find that ε * ≈ 7.1. This value of the intermediate dielectric constant is reasonable; it is close to the dielectric constant at the frequency separating librational and orientational fluctuation modes in water at ∼5, but of course this value is generally poorly known.

III.A. Indirect Solvent Contribution to the Hydration Force.
In the spirit of Marcělja-Radic, 9 we can calculate the indirect solvent contribution to the hydration force by finding the derivative of the energy obtained from our Hamiltonian. We thus consider a film of water bound by two planar interfaces, infinite in the xy plane, located at z = 1 and z = L − 2 . Rather than explicitly placing the planes at z = 0 and z = L, i.e., 1 = 2 = 0, we define the position of the surface (where the boundary conditions are placed) as a variable, which will let us take roughness into account later on.
The symmetry of the problem allows us to consider only the longitudinal polarization depending on z, such that P(r) = P(z) u z . This simplifies the functional to one dimension, thus yielding the energy per unit area. Calculating the derivative of this 1D functional will give us the hydration pressure. We obtain the longitudinal polarization density profiles in space by variational minimization of the Hamiltonian, resulting in a pair of coupled differential equations in P 1 (z) and P 2 (z), which are solvable with the following boundary conditions in P 1 (z) and P 2 (z): These boundary conditions denote a system where water is polarized at two similar surfaces, indicated by the opposite sign of the total polarization density, P̅ 0 at each surface. However, the contribution to P̅ 0 from each polarization mode is not immediately clear. Hence, in the boundary conditions, we introduce the parameter α, which divides the total polarization density at the surface between the two modes, such that The Journal of Physical Chemistry C pubs.acs.org/JPCC Feature Article bound charge at the surface, ρ ̅ b , should be. Hence, for all calculations performed, the value of ρ ̅ b taken is that which minimizes the energy (eqs 1−3). Solving numerically the coupled equations in eq 15 with eq 16 as boundary conditions, we plot the polarization density profiles between two uncharged surfaces in Figure 2. It is clear that the decay and oscillation in the profiles are determined by the correlation lengths, the values of which are obtained from fitting the dielectric response functions. The resonance peak in χ(k) manifests in P 1 as decaying oscillations with periodicity λ o (corresponding to the diameter of a water molecule) and decay length λ d, 1 . Similarly, the Lorentzian behavior in the small wavenumber region of χ(k) manifests as exponential decay of P 2 with characteristic length λ d,2 . It is not purely exponential however�some waving in the profile is present due to the coupling between P 1 and P 2 (cf. eq 15).
Injecting the numerically calculated P 1 (z) and P 2 (z) profiles into the Hamiltonian (eqs 1−3) and performing the integral over z, we obtain the energy per unit area as a function of the surface separation from which we can calculate the solvent contribution to the force (per unit area), Π(L), between the two surfaces. This is given by the derivative of the energy per unit area w.r.t. the distance between the surfaces, Obtaining the force as a function of the positions of the boundary allows us to introduce the overall effect of surface roughness. We do this by averaging the hydration force over distributions representing the two surfaces: where ρ( ) describes the surface distribution. Here we use a truncated Gaussian distribution for each surface, normalized over half the separation distance, L: 19) where is the mean position of the surface, δ is the half-width of the distribution and Θ is the Heaviside step function. It should be noted here that any arbitrary surface distribution, e.g., a uniform distribution, yields qualitatively similar results.
The computed hydration force profile and the effect of smearing the force is shown in Figure 3. For sharp interfaces, we see an oscillatory force profile, but when the surface is smeared, the oscillations rapidly disappear. A surface roughness corresponding to the size of a water molecule (2δ ≈ 2.5 Å) is enough to completely remove oscillations from the profile, leaving only an exponentially decaying force. It is clear from comparing this force with the one obtained in the Marcělja-Radic model (see Appendix B) that it is not just the Lorentzian mode that survives after smearing. There is some additional contribution to the remaining force, likely a result of the coupling between the exponential correlations present in the Lorentzian and overscreening modes, the latter disappearing upon smearing in the uncoupled case.
Interestingly, in the case of purely exponential correlations (the only 'survivors' after surface smearing), when the boundary conditions are such that the polarisation density vector at each surface is oriented in the same direction along the z-axis (i.e., have the same sign), we see attraction rather than repulsion. This case may have its own interesting features 19 that require further investigation.
III.B. Interaction between Mica Surfaces. The above theory treats only the correlations we observe in pure water polarized by a surface. In reality, charged surfaces and electrolyte subsystems will couple with the solvent structure, and so we cannot strictly reproduce experimental data with a superposition of direct (DLVO) and indirect contributions. However, in the classical experiments of Pashley and Israelachvili, 25,27 the ion  Figure 1. The resonance peak in the dielectric response function indicative of overscreening and oscillations can be seen present in P 1 (z), whereas P 2 (z) is a result of Lorentzian behavior of the dielectric function at low wavenumbers. Importantly, the profiles decay to bulk properties over the distance of 1 nm. concentrations are small, such that the Debye length is much greater than the correlation lengths in pure water. For these systems, decoupling the solvent and electrolyte contributions to the force is a decent approximation to make. The direct local electrostatic force between two flat surfaces at large distances in aqueous solution has previously been well described by DLVO theory, and the solvent contributions at short distances can be included using the theory developed in this paper. Hence, to reproduce experiment, we can use the following approximation: Often in crossed-cylinder SFA experiments, the force is normalized to the cylinder radius (F/R). This can be shown to be equivalent to 2 using the Derjaguin approximation, is the interaction energy per surface area between the two cylinders. 8,50 Hence, the DLVO contribution here will be given by the classic formula for the repulsive free energy per unit area between two planar surfaces, 51 where n b is the bulk electrolyte ion density, κ D = λ D −1 is the inverse Debye length, q i is the valence of the electrolyte involved (q i = 1 for monovalent electrolytes), q e = 1.6 × 10 −19 C is the elementary charge, and ϕ 0 is the electrostatic potential at the surface. Given the maximum polarization density of water P̅ 0 = P̅ 1 + P̅ 2 ≈ 30 μ C/cm 2 at the interface, we assume a significant fraction of water is oriented in the same direction at the surface; hence, we use P̅ 0 = 10 μ C/cm 2 below. With this we faithfully reproduce the profile found by Pashley and Israelachvili (see Figure 4). The magnitude, decay lengths, and characteristic lengths seem to coincide very well, speaking in favor of the above developed picture. When we smear the solvent contribution (δ = 0.15 nm), we qualitatively obtain the profile found from earlier experiments between mica surfaces. We see that the decay length of the smeared force (λ d,2 = 3.65 Å) does not exactly coincide with the experimental value (λ d = 8.6 Å), 27 which can be for a number of reasons. This decay length strongly depends on the fit of the minimum in the response function. Higher accuracy and resolution in this region would result in a more informed fit and hence more accurate correlation lengths.
Additionally, we have neglected any coupling between the solvent correlations and lateral correlations along the surface. Including these would account for variations in the decay length depending on surface properties, as observed for lipid bilayers. It is also possible that these lateral surface correlations could smear the water structure and hence the hydration force without the need to include this surface roughness.
However, what we unambiguously show here is the necessity of the Lorentzian component to the hydration force to retain an exponentially decaying force profile between rough surfaces. Oscillations are introduced by including the overscreening resonance in the response function, and they survive only when the surfaces are smooth and sharp. With rough surfaces, these oscillations end up being suppressed due to dysphasia.
Having focused on the role of polarization fluctuations, one must not forget the importance of the fluctuations of molecular density�the focal point of a set of works of David Chandler 52,53 and Daniel Borgis. 54 In the context of our study, it is important to note that the resonance peak in the response function emerges naturally as a result of coupling between the polarization and density fluctuations; 37 see also ref 55 (where, however, the ability to reproduce the overscreening peak was imprinted already in the polarization part of the free energy functional). Using more complicated free energy functionals that involve both polarization and density variables may help to describe systems where each of those variables may be independently perturbed through the boundary conditions, which may be needed in the interpretation of a number of phenomena, including hydration forces. 56 Doing this here would have however complicated the analysis, and for the sake of transparency, we leave it for future work.
As we have not explicitly considered the role of density correlations here, we should, strictly speaking, refrain from conclusions about the solvation forces mediated by nonpolar liquids. Indeed, it is well-known that the force profiles between surfaces separated by nonpolar media also exhibit oscillations, as seen in experiments, 57−59 as well as in theoretical hard sphere liquid studies. 60,61 However, it has also been established in liquid/vapor interface theoretical studies (also with hard sphere models) that the capillary-wave-induced "roughening" of the interface suppresses 62 or smears 63 the oscillations in the hard sphere density profiles, when averaged along the Gibbs dividing surface. So far, we have considered the system as a decoupled superposition of the double layer and solvent structural contribution, which is only valid for very low electrolyte concentrations. However, we must be able to couple the two structures together and see whether this correlated structuring of water due to the interface does affect the structure of the double layer and how smearing influences this structure. An attempt to do this is made below in the context of an electrical double layer which, as we will see, leads to remarkable results.

IV. DOUBLE LAYER STRUCTURE AT SMALL ELECTROLYTE CONCENTRATIONS
After this study of hydration forces and water polarization profiles confined between two surfaces, it is tempting to capitalize on the obtained findings and investigate how the structure of the electrical double layer is affected by these features. Independent of the intrinsic logic of this study, it has strong motivation. There is currently a lot of debate as to how water structuring could influence ion distribution near the interface. The interest to this question however has long history, based either on the earlier exploration of nonlocal electrostatic effects described within the concept of Lorentzian correlations 64 or on more general formalism 65 which describes how this structuring of water will influence the way ions are distributed close to the surface. There already exists a significant body of work, the most notable of which includes the Gouy−Chapman and Poisson− Fermi approaches, which will be expanded upon in the following sections when relevant. Here we develop a Hamiltonian for the double layer which couples the potential, polarization density, and ion density profiles.

IV.A. Poisson−Fermi Free Energy Functional.
To continue in the vein of a continuum theory, we extend the widely used Poisson−Boltzmann (PB) approach to describe the electrolyte in this system, coupling it to the water structure. PB theory is very successful in predicting ionic properties near planar and curved surfaces, and the resultant double layer forces. However, it does not include finite size of adsorbing ions, i.e., the excluded volume effect, and so it can drastically overestimate the ion concentrations close to the surface. This overestimation of ion density could cause problems when incorporating oscillating polarization density profiles into the model. Such strong oscillations could draw ions into the deep potential wells producing very large, physically inaccurate ion densities. It is therefore sensible for us to adopt a modified PB approach which accounts for these steric effects. This approach is often referred to as Poisson−Fermi due to the Fermi−Dirac-like nature of the ion distributions, and it has been implemented a number of times in the past for electrolyte/ionic liquid systems. 66−69 Within the mean-field approximation, we can write the total free energy of the system as a functional in terms of the electrostatic potential ϕ(r), the ion densities n ± (r), and the two orientational polarization density modes P 1,2 (r) as described above. This free energy functional can be broken down into electrostatic, entropic, and correlation terms: At first glance, it may appear that there is an inconsistency between this electrostatic free energy functional and the one introduced in eq 2 for pure water, in that eq 2 has omitted relevant terms in ∇·D (see ref 70). However, in the case of pure water, D = 0 or is constant for charged plates, and so all terms in ∇·D are vanishing, and both models remain consistent with each other.
The correlation term is approached in the same way as above in eq 3, as a Landau−Ginzburg expansion of two orientational polarization modes. We describe the entropic contribution to the free energy by a lattice gas approach, where the solution is modeled as a cubic lattice with each cell of volume a 3 filled with one ion or water molecule, such that a 3 (n + + n − + n w ) = 1. For simplicity, a is taken to be the same molecular size of all ionic species and solvent. Expressing n w in terms of n + and n − , we obtain k T n n a n n a a n a n a n a n a where k B T is the thermal energy. Here the entropies of the positive (n + ) and negative (n − ) ions are stated in the first two terms, and the last term explicitly states the entropy of the water molecules (n w ). To find the potential and polarization profiles between the two surfaces, we first approach the ion density profiles by variational minimization of the functional with respect to n ± , and coupling the system to a bulk reservoir, such that n / tot = ± ± , where μ ± are the chemical potentials of the ions. In the bulk, for a 1:1 electrolyte, n + = n − = n b , where n b is the bulk electrolyte ion density. We define the bulk volume filling fraction, ξ b = n b a 3 , and hence we can write an expression for μ ± , This yields the ion density profiles: where the potential and the ion densities are scaled by Φ = q e ϕ/ k B T and n ± * = n ± /P̅ 0 , and hence a* 3 = P̅ 0 a 3 . Following the scaling, such that P i * = P i /P̅ 0 , variational minimization and In the absence of ions, i.e., ξ b = 0, we recover the system for pure water in eq 15�see Appendix C for specific details on obtaining potential profiles in pure water. Similarly, if we allow the coefficients of the higher order gradient terms in the polarization density expansion to equal zero such that the solvent behaves with a bulk dielectric constant, and in the limit of ξ b ≪ 1, the Poisson−Boltzmann equation for a 1:1 monovalent electrolyte is recovered, ∇ 2 Φ = κ D 2 sinh Φ. As above, the symmetry of the problem allows us to take this in just the z-direction, greatly simplifying the higher order operators. To solve these equations for the potential and polarization density profiles, we must carefully consider the boundary conditions. To recreate a system with chemisorbed water at the surface, the boundary conditions for polarization density remain the same as above, see eq 16. A more sophisticated way to approach the polarization density boundary conditions could be to introduce an induced P̅ 0,ind alongside a chemisorbed P̅ 0,chem polarization density to allow variations in the water structure at the surface according to the surface charge density. However, if the chemisorption is reasonably strong, we can assume that the contribution via the induced part is small (valid for smaller surface charge densities), and we can neglect this here, simplifying the problem. For the electrostatic potential, we use Neumann boundary conditions, which are determined here by the surface charge density and the polarization density at the surface. From the definition of the displacement field in electrostatics, D = ε 0 E + P, we see that at the surface, D = σ ̅ 0 and P = P̅ 0 , and using E = −∇ϕ, we obtain for the boundary conditions at the surface: It is important to note that this formulation is designed to treat only small concentrations of electrolyte. Linearizing the system of equations in (28) and calculating the response functions, we obtain for the dielectric function here, ε c (k): where ε(k) is the dielectric function of pure water as described in eq 10. It is clear that increasing κ D , can violate the DKM rule for ε c (k). To be able to describe the domain of higher electrolyte concentrations, a more sophisticated description of the ionic subsystem would be needed.

IV.B. Water-Mediated Double Layer.
We reiterate that the attention to the interplay of the molecular nature of the solvent and the statistical mechanics of the ion distribution in the electrical double layer is experiencing a renaissance nowadays, 71,72 but the problem per se is not new. We refer the reader to a detailed review of related effects in the monographic chapters in refs 73 and 74, in particular to the section 2.1.10 devoted to this issue. Already then it was clear based on the results of molecular simulations that simple pictures of the ionic distributions near the charged electrodes are very much different from the standard double layer theories of the Gouy−Chapman type and its primitive solvent model extensions. Below, we explore this issue in full detail based on the formalism presented above, staying away from the case of concentrated electrolytes which would inevitably require application of a more involved description of the ionic subsystem. We start first with the findings displayed in Figures 5 and 6, obtained respectively for ideally sharp and for slightly smeared surfaces: (i) Even in a very diluted electrolyte, ions do get trapped in the potential wells created by the overscreening dielectric response of water to the charge of the electrode, known also as layering of water molecules, as we see in the first row of panels in Figure 5. Such "water wells" are deep enough to keep ions in, against the entropic drive to spread them around. (ii) Ions do not follow the potential profiles due to pure water literally, as their presence disturbs their corresponding potential wells, albeit not much. Ion concentration in the wells is limited by the effect of "excluded volume" (finite ion size), so that in the first 3−4 wells you actually have the same, maximum possible concentrations of cations or anions, as seen by the plateauing profiles in Figure 5a. We thus cannot see the effect of concentration of ions in the bulk, as the occupation of the wells is very much close to the limit determined by the excluded volume of ions; we may call it a Langmuir-type electrosorption. (iii) However, as we see in Figure 6, a small smearing of the interface, which as above we saw suppresses the overscreening resonance effects in pure water, eliminates the ion density oscillation peaks, leading to an amazing result which looks like a Stern layer near the surface alongside a classical Gouy−Chapman layer. This may be a reason, in many cases, why the good old Gouy− Chapman-Stern theory works so well! (iv) Fixed boundary conditions on the water polarization at the surface not only shift the potential of zero charge, but affect the overall distribution of ions. Indeed, larger values of polarization density at the boundary, which cause more persistent overscreening, stimulate adsorption of ions into the "water wells", so the oscillations in the ion distribution become more pronounced at longer distances from the interface. (v) The effect of the absolute value of electrode potential is profound: it not only Gouy-like compresses the electrical double layer but also suppresses the overscreening. The series of plots shown in Figures 5 and 6 demonstrate these findings. In these, we examine the ion and electrostatic potential distributions from a single positively charged surface. In all cases, the distributions are compared against the Gouy−Chapman result, where the surface potential is calculated via the Grahame equation. All mathematical details for Gouy−Chapman theory can be found in Appendix D. Of course, polarized pure water near a charged interface establishing oscillating electrostatic potential profiles and oscillating profiles of counterions and co-ions at electrodes are not new findings. Such double layer structures have been seen in numerous molecular dynamics and classical density functional simulation papers; see, e.g., refs 75−80. However, most of these papers considered concentrated solutions or ionic liquids. As such, they have not generally established the precise relation between the water electrostatic potential profile and the density profiles of cations and anions. Two interesting recent papers 71,72 have come closest to this point in considering moderate electrolyte concentrations (ref 72 unfortunately did not specify the exact value of salt concentration in their simulations); the reported profiles in both papers seem to support our findings.
There have been experimental studies pertinent to the subject of this article, performed not with SFA but with AFM. Of particular interest to us is the work of the Garcia group, 81 which uses the 3D-AFM experimental techniques developed by the Fukuma group. 82 Garcia's 3D-AFM results at the spontaneously negatively charged mica surface in contact with electrolyte interface, unlike some other works along these lines, are not purely restricted to the case of highly concentrated salt solutions. That group also presents results down to 10 mM concentration. Even at such a low concentration, they show evidence of distinct layers of water (hydration layers) and cations in the normal direction to the surface. This is in line with the results we have obtained for a sharp interface ( Figure 5). This correspondence is remarkable, but still should be taken with a pinch of salt as alluded to in the work by the De Yoreo group, 83 because the AFM tip may in principle substantially perturb the native interfacial structure of solution. However, if the probe tip is smaller than the characteristic length scale of roughness of the surface oscillating structure in the liquid has been observed. 84 Surprisingly, this has been reported even when measuring forces against materials like silica, 85 infamous for their roughness. 86 Interestingly, all the above-mentioned works also explore the role of lateral inhomogeneity of the sample surfaces. While the latter is not considered in our paper, previous work that studied the interplay between the correlations in the XY plane with Lorentzian dielectric response in water shows that the characteristic correlation length of the surface couples with the exponential decay length of water leading to a new, reformed characteristic decay length. 17 While that theoretical study was applied to lipid bilayers, the results can be extended along these lines onto the case of solid surfaces. At present it is not clear how lateral correlations would couple with the overscreening effects, and this would be an interesting avenue for further work.
Although there are many new experimental methods that could resolve the structure of the interface with atomic resolution, 87,88 as well as methods of molecular simulations, there are not enough systematic studies of the distribution of cations and anions in electrolytes of low concentration at different electrodes and at different electrode potentials. Work in this direction started at the Argonne National Laboratory years ago. 89 Recently, important results in this direction have been obtained by Lee et al., 90 in which using X-ray scattering, the salinity-dependent EDL structure at negatively charged mica− water interfaces was studied in the interval of concentrations from 0.022 to 6.7 M, also supported by the corresponding molecular dynamic simulations. The effects of overscreening were observed, the ion layering attributed to strong ion correlation effects. Obviously at high concentrations it will be difficult, if not impossible, to discern the effects of correlations in water and in ion "clouds", as the two subsystems are intertwined. As ref 90 did not, however, focus on correlations between the water structure and the structure of double layer in the most diluted solutions that they have studied, we cannot directly compare those aspects of our findings with them, but the observation of decaying ion oscillation profiles near the interface in diluted solutions is what our field-theoretical model also suggests. The authors of ref 90 note that the EDL structure, characterized by the development of both lateral positional correlations between adsorbed ions and vertical layering of alternating cations and anions, is reminiscent of the structures of strongly correlated ionic liquids at electrified interfaces. The iconic data for the ionic concentration profiles of ref 91 obtained on a spontaneously charged smooth sapphire surface does show similar overscreening and layering, but the periodicity of the latter is, of course, different from the one observed in ref 90, where water and hydration plays a role in the structuring of ionic distributions.
Is the trapping of ions in diluted solution in water "wells" a real phenomenon, or is it exaggerated in our model? If it exists, does this effect weaken for intermediate electrolyte concentrations, as a result of some intereference between ion and water correlations? At present, we can only speculate about it. There are two at play here: ionic correlations and water correlations. Indeed, there should be a coupling between the two, but each effect will be dominant at different concentration limits. At low concentrations, the correlated water structure seems to play the primary role in the oscillatory structure of the double layer (for sharp interfaces), although the immediate environment of each ions may be perturbed by their presence. 35,92 At high concentrations, ion correlations (not described in our model) take charge, also leading to oscillatory behavior, seen even in the primitive solvent description. 93 However, for intermediate concentrations, it is not exactly clear how these effects would interfere/couple with each other. If the two effects do quench each other, at how large or small an ion concentration would this quenching occur?
It should made clear that in our model we do not consider any specific ion adsorption at the surface, usually accounted for by special local terms in the free energy. Hence, the only parameter involved in controlling the ions is the bulk ion concentration, and it is this value that controls the accommodation of ions in the "water wells". When the ion gets trapped in the first "water well" near the interface, it, however, looks as if it is adsorbed with its hydration shell at the interface. By classical definition in electrochemistry, such ions are called "nonspecifically adsorbing" and are not chemisorbing because they do not take part in any electron exchange with the electrode. Upon formation of that first layer of hydrated ions, at their dense packing, the average distance between them would be of the order of 0.8 nm. The corresponding value for the ion coverage should stand for the maximum concentration (saturation value) in the first layer of hydrated ions, and such values should be considered in the expressions for excluded volume term in the free energy.
IV.C. Double Layer Capacitance. While the plots above show predictions for how water structure contributes to the double layer, they are difficult to validate as few experiments have been performed looking at ion profiles in such dilute systems. Despite our model receiving some validation from our study into hydration forces, the interplay between surface charge, water, and ions can also manifest in another important measurable interfacial quantity in electrochemistry, the double layer capacitance, giving our model another route to explore. The double layer capacitance integrates all the above-discussed effects. As we shall see, the theory gives results for the capacitance comparable with experiments without the need to invoke any compact (Stern) layers. A great deal of effort in the past has been made to rationalize results obtained for the capacitance, and the consequent picture they lead to for the double layer. 73 Now we will explore what our description suggests for it.
The total differential capacitance of the double layer, C, is defined as where ϕ 0 is the surface potential. Using D = ε 0 E + P at the surface, Here we assume that P̅ 0 is independent of Φ 0 . Such an approximation is usually associated with strong chemisorption of the water at the surface, but we will keep it here for a general case as a matter of simplification. We therefore simply obtain the standard result for the capacitance: It should be noted however that this does not remove any effect of polarization density on the capacitance. As the potential is calculated from the coupled equations in eq 28, there is an inherent dependence on the polarization density. Using eq 33, we plot the differential capacitance as a function of surface potential for different bulk concentrations n b and surface polarization densities P̅ 0 . For P̅ 0 ≈ 0 μC/cm 2 , we observe a well-defined symmetric capacitance minimum. However, instead of growing exponentially at large potentials as in Gouy−Chapman theory, it decreases at high potentials at the "wings". This is similar to the established Poisson−Fermi approach which takes excluded volume into account (for details and relevant mathematical expressions plotted in Figure 7, see ref 69). However, we should make a note that the behavior in our model at these higher potentials may be inaccurate as a result of the approximation that P̅ 0 is independent of ϕ 0 . Despite this, at smaller electrode potentials where the model is more likely to be valid, the capacitance shows a broader minimum, which is more in line with experimental results. 94 As with standard Gouy− Chapman theory, we also observe an overall increase in the capacitance with concentration, however we will address the differences in the model regarding this in the next section on Parson−Zobel plots.
An interesting observation to note is the effect of P̅ 0 . In particular, it introduces asymmetry into the system, and a shift in the value of the potential of zero charge (pzc). In experiment, this coincides with the potential at the Gouy−Chapman minimum. We cannot directly compare values of the pzc with experiment, as the potential drop here is measured relative to the bulk as opposed to a reference electrode, but we show here that the strength of water chemisorption to the metal surface plays a significant role in the positioning of the capacitance minimum. Additionally, a nonzero P̅ 0 introduces further features such as shoulders and secondary minima which strengthen with concentration, as observed in experiment. 95 where within Gouy−Chapman theory, C GC is given by where C D = ε 0 ε b κ D is the Debye capacitance. We have shown above in the study of the ion layer profiles that this Stern layer does indeed arise from structured water hydrating the surface without the need to explicitly invoke one. This result is interesting on its own, but we can go one step further by performing a Parsons−Zobel (PZ) analysis on the model. Here, the inverse capacitance of the model at varying concentrations is plotted against the corresponding inverse capacitance predictions of Gouy−Chapman theory. It is clear from eq 34 that a plot of C −1 against C GC −1 yields a linear dependence with gradient equal to 1 and intercept equal to C H −1 . Many experimental studies with metal electrodes, Hg in particular, have been very faithful to this, serving almost as a verification of the GCS capacitance model 96 and, hence, the simplified picture of the double layer. However, more recent studies with Au(111) and Pt(111) have shown significant deviations to this behavior, indicating a much different picture to the Gouy−Chapman one. Ojha et al. observed that for these surfaces, a PZ plot does yield a linear relationship, but with a drastically reduced gradient of 0.179 for Au(111), and 0.041 for Pt(111). 97 It was suggested by Parsons and Zobel that this observation may arise from specific adsorption of ions at the surface, leading to a higher measured capacitance than predicted. This motivated a theoretical treatment by Schmickler where this effect was explicitly taken into account, 98 resulting in a reduced PZ gradient in correspondence with experimental results for Ag(111) electrodes. 99 Further theoretical work has been done since, introducing an effective spatially dependent dielectric constant near the electrode. 100 We also approach this problem of the PZ plot gradients, however from a slightly different angle. We see from the ion profiles above in Figures 5 that structured water near the sharp interface results in accumulation of ions at the interface, akin to adsorption of hydrated ions. This begs the question, could the chemisorption and hence preferential orientation of water at the surface give rise to this deviation from ideal GCS behavior?
In Figure 8, we plot the inverse capacitance in our model at the pzc for a number of ion concentrations against the corresponding inverse Gouy−Chapman capacitance, C GC −1 . It is not required to use the more general Poisson−Fermi capacitance here, as we measure the capacitance at the pzc, where the two models are essentially equivalent. A linear relationship is observed with deviations at higher electrolyte concentrations. This is not a new finding, and has previously been attributed to solvent structural effects. 101 Perhaps the more striking feature in this plot is the effect of P̅ 0 �a small increase in the polarization density of water at the surface leads to a significant decrease in the gradient, in line with experimental results found by Ojha et al. 97

V. CONCLUSIONS
Using a consistent field theory model of water and aqueous solutions between electrified interfaces (charged and polarized) we have demonstrated that hydration forces can contain oscillating features depending on whether the surfaces are ideally sharp or smeared. Oscillations in forces are caused by water layering and overscreening. The predictions of the theory are consistent with the known experimental observations. The (successful) attempt to describe hydration forces for a slightly smeared surface has served as indirect evidence of the Lorentzian contribution to the nonlocal dielectric function, the existence of which is still a subject of debate (some accepting it by default, some neglecting it altogether).
We have shown that the same effects cause oscillating profiles of cations and anions near the interface even in a very diluted solution. Again smearing of the interface almost fully removes these oscillations after a couple of water layers, and the ion distribution and the potential distribution is well reproduced by the classical Gouy−Chapman theory with a correction�a contribution of a compact layer kind, but determined by the polarization of the first few water layers. Although this may be received with relief by many electrochemists, still the fine structure of layers nearest to the surface is a question of first rate importance, for instance for electrochemical kinetics.
We have shown various signatures of these fine water structure effects in this paper coming from the theory, most of which will await verification by scattering and AFM experiments. Our present article is intended to further stimulate such experiments.
We also present the results for capacitance of a semi-infinite double layer and compare them with classical, as well as recent experimental results. The analysis unravels the contribution of the nearest layers and the interface into this fundamental integral characteristic of electrochemical interfaces.
Most importantly for many, the tiny smearing of the interface results in a picture looking formally similar to the Gouy− Chapman−Stern one, when the first peak of counterions sits near the surface at a distance of one water-molecule diameter, and there are no other peaks, but a standard Gouy−Chapman decay of the charge density distribution down to the electroneutral bulk.
A final, but important note: the specific nature of ions is not included in the considered model. First, the description presented above is based on linear response approximation (apart from the description of the double layer), and so nonlinear effects such as "positive" or "negative" hydration of ions (of kosmotropes or chaotropes), 35,92 required for description of different alkali halide ions have not been accounted for. In particular, the oscillatory component of the hydration force has been observed to be ion specific, being present for strongly hydrated (kosmotropic) ions such as Li + and Na + , whereas they are completely suppressed in the presence of weakly hydrated (chaotropic) ions such as Rb + and Cs + . 102 Introducing this ion specificity in this principal study would complicate, and potentially obscure, the picture, as our desire was to describe a reference system of effects beyond nonlinear ones. As we know that the pictures based on the nonlocal dielectric response of the bulk water may be substantially distorted by the mere presence, charge, and size of ions, 35,92 we cannot be fully certain that this "reference" system will be instructive enough. However, analysis of the model in this paper shows that it at least captures some of the basic effects, although imposition of the preferences of water structuring the ionic distributions may be partially overestimated. It would however be interesting to use the model we have formulated as a reference result with which we can compare against these nonlinear effects.

A.1. Molecular Dynamics Simulations
We simulate a cubic water box of side size L = 6.5 nm composed of N w = 9033 water molecules.
Simulations are performed using GROMACS 2021 molecular dynamics simulation package, 103 where the integration time step is set to Δt = 2 fs. Simulation boxes are periodically replicated in all directions, and long range electrostatics are handled using the smooth particle mesh Ewald (SPME) technique. Lennard-Jones interactions are cut off at a distance r cut = 0.9 nm. A potential shift is used at the cut-off distance. All systems are coupled to a heat bath at 300 K using a v-rescale thermostat with a time constant of 0.5 ps. We use MDAnalysis to treat the trajectories. After creating the simulation box, we perform a first energy minimization. We equilibrate the system in the NVT ensemble for 200 ps, and afterwards in the NPT ensemble for another 200 ps using a Berendsen barostat at 1 bar. Production runs are performed in the NVT ensemble for 20 ns.
We performed simulations with SPC/E water model, a threepoint charge model. In this case, the oxygen atom carries both LJ center and the negative charge.

A.2. Longitudinal Susceptibility
To compute the k-dependent susceptibilities, we use the fluctuation-dissipation theorem, relating the response functions to the polarization fluctuations as follows: Using the relation ρ b (r) = −∇·P(r), one can express the longitudinal susceptibility as a function of the charge structure factor S(k) and gets The charge structure factor in the Fourier space can be decomposed into an intramolecular and intermolecular part, 22 where d IJ is the intramolecular distance between atom I and J. At low k the precision of this expression of the structure factor become pretty low as the function h(r) is obtained on a finite range imposed by the box size. To solve this problem, we proceed as follows. For k < 2.5 Å −1 , we take into account the periodicity of the system, calculate the charge structure factor for d i s c r e t i z e d v a l u e s o f t h e w a v e l e n g t h k , k L n n n 2 / where k is a vector and d AiAj stands for r⃗ Ai − r⃗ Aj .

A.3. Statistical Treatment
For the longitudinal susceptibility at low k, the error bars are derived following the reblocking method. 104 ■ APPENDIX B. SMEARED MARČELJA-RADICH

YDRATION FORCES
Here we derive the necessary expressions required to calculate the hydration forces between two smeared interfaces in the Marcělja-Radic̀(Lorentzian) model. Using our {K ij } parameter formalism, we started with a simple quadratic LG expansion in the polarization density to the first gradient term: We see that following the method in the main text, we obtain for the dielectric function: Hence, we see that We also see that, from the pole on the imaginary axis ε MR (k), that we have a correlation length, K K / 22 21 = . In comparing the two models in Figure 3, we use ξ = λ d,2 , the fitted pure exponential decay length from the hybrid model presented above in the main text.
Considering the symmetry of the problem as above, we only need to consider the normal (z) direction, as P(r) = P(z)u z . Variational minimization gives the following differential equation to solve for P(z): We solve this between two surfaces placed at z = 1 and z = L − 2 , where L, 1 and 2 have the same meaning as in Section. III.A. Fixing the polarization density in opposite directions at each surface, such that P(z = 1 ) = −P(z = L − 2 ) = P̅ 0 , and solving, we obtain the polarization density profile, P(z). In the same way as above, we reinsert P(z) into eq B1, and calculate the force using We can then smear this using eq 19, to have direct comparison with the smeared forces in the hybrid model.

PROFILES BETWEEN OPPOSITELY CHARGED PLATES IN PURE WATER
Despite having formulated a theory for pure water resulting in polarization density profiles in Sections II and III of this paper, it is not immediately clear how we could obtain electrostatic potential profiles. Indeed, this can only be done between equally oppositely charged surfaces to preserve electroneutrality, and so with the formalism established in Section IV we can derive the electrostatic potential in the gap. As stated above, setting ξ b = 0 will reduce the system of equations in 28 to the equations for pure water in 15. To calculate the electrostatic potential, however, we need to consider the line: Integrating twice, and using the boundary condition as required in eq 29, we find the potential profiles due to pure water q P k T z P z P z P z c d ( ) ( ) from the requirement that Φ(z = L/2) = 0 in between the two plates. Hence, solving symbolically the coupled differential equations in eq 15 in the main text allows calculation of the potential distribution from the pure water polarization density profiles from eqs B2 and B3. = is the inverse Debye length. Solving this separable differential equation with the boundary conditions ϕ(z = 0) = ϕ 0 , where ϕ 0 is the electrostatic potential at the surface, after a small amount of manipulation, we obtain the Gouy−Chapman potential distribution where Φ 0 = q e ϕ 0 /k B T. Using the relation from bulk electrostatics, dϕ/dz| z=0 = −σ ̅ 0 /ε 0 ε b , where σ ̅ 0 is the surface charge density, we simply obtain the well-known Grahame equation from eq C2: ■ REFERENCES