Water molecules mute the dependence of the double-layer potential pro ﬁ le on ionic strength

We present the results of molecular dynamics simulations of a nanoscale electrochemical cell. The simulations include an aqueous electrolyte solution with varying ionic strength ( i.e. , concentrations ranging from 0 – 4 M) between a pair of metallic electrodes held at constant potential di ﬀ erence. We analyze these simulations by computing the electrostatic potential pro ﬁ le of the electric double-layer region and ﬁ nd it to be nearly independent of ionic concentration, in stark contrast to the predictions of standard continuum-based theories. We attribute this lack of concentration dependence to the molecular in ﬂ uences of water molecules at the electrode – solution interface. These in ﬂ uences include the molecular manifestation of water ’ s dielectric response, which tends to drown out the comparatively weak screening requirement of the ions. To support our analysis, we decompose water ’ s interfacial response into three primary contributions: molecular layering, intrinsic (zero-ﬁ eld) orientational polarization, and the dipolar dielectric response.


Introduction
][7][8][9][10][11] Certain features of the interfacial potential prole can be systematically modied, for example by changing the applied potential or the chemical composition of the electrolyte solution. 5,7,10,12In this manuscript, we study the effects of changing electrolyte concentration on the interfacial potential proles, specically comparing the predictions of continuum-level theory and all-atom molecular dynamics (MD) simulation.We highlight that the potential prole derived from MD simulation exhibits little apparent dependence on ionic strength, in stark contrast to the predictions of the standard textbook theory of Gouy and Chapman.We demonstrate that the source of this discrepancy is that the continuum theory does not account for the molecular effects of the solvent and its polarization response.Our simulation data reveals that these effects play a dominant role in shaping the interfacial potential prole.By analyzing the simulations of neat liquid water between charged electrodes, we attempt to deconstruct water's various contributions to the interfacial potential prole.These contributions include oscillations from molecular layering, a symmetric potential drop between each electrode and the bulk from the inherent orientational bias of water-metal interactions, and a modied dielectric response due to excluded volume effects at the electrode boundary.
It has long been appreciated that the performance of electrochemical systems can be systematically modied by the concentration of ions in the supporting electrolyte.Since catalytic electrochemical reactions oen proceed through surface intermediate species that reside in the electrochemical double layer (EDL), tuning screening behavior in the EDL (e.g. by changing the ionic strength of the solution) can meaningfully change electrochemical reaction rates and mechanisms.In the case of reactions with outer-sphere mechanisms described by Marcus kinetics, theoretical, 13,14 and experimental, 15 work has demonstrated that the reorganization energy (and hence reaction rate) is sensitive to the screening length-scale set by the electrolyte ionic strength.For inner-sphere electrochemical mechanisms, the effect of changing the ionic strength on the reaction rate can be more difficult to model, yet there are a number of experimental studies that have experimentally characterized the effects of changes in electrolyte composition and ionic strength on hydrogen evolution and CO 2 reduction kinetics. 16,17stablishing the boundaries of validity for standard theories and models, and identifying the molecular origins of their failures, is vital to our ability to advance the eld of electrochemistry.Here, we utilize MD simulation to explore water's role in mediating the screening response of the aqueous interface.We compare our results to the predictions of standard theories, which generally treat liquid water as a dielectric continuum.Since molecular dynamics simulations fully represent the molecular structure, e.g., size, shape, and orientational correlation, they are well suited for validating the assumptions that underlie common theories and models.The manuscript is organized as follows.In the following section, we review double-layer theory.Then, in Section 3, we present results exploring the effects of changing ionic strength on the electrostatic potential prole of the interface.In Section 4 we evaluate the specic roles that water plays in mediating these effects and in shaping the potential prole more generally.Finally, following a brief conclusion, we describe our theoretical and simulation methodology.

Theoretical descriptions of the electric double-layer
8][9] Theoretical descriptions of EDL structure and its role in screening phenomena have a rich history.9][20] Under this assumption, as illustrated in Fig. 1A, the EDL is a microscopic parallel-plate capacitor, with a constant electric eld between the electrode and the OHP, and no variation in the potential from the OHP onwards.This EDL theory of Helmholtz and Perrin offers important physical insight into the role of ions in screening elds originating from electrodes, but has signicant limitations.
The complete screening layer of Helmholtz and Perrin is entropically unfavorable and therefore not stable under standard thermal conditions.Entropic effects cause the screening layer to spread, broadening the width of the EDL.Physically reasonable electrostatic screening congurations therefore feature a diffuse cloud of neutralizing ionic density.An elementary description of potential variations in the EDL that accounts for such entropic effects originates from the Poisson-Boltzmann equation, 21,22 where 4 is the electrostatic potential, c 0 is the bulk ion concentration, e is the fundamental charge, 1/b = k B T is the Boltzmann constant times temperature, and 3 and 3 0 are the continuum dielectric permittivity and the vacuum permittivity, respectively.The summation runs over all ionic components of the solution, indexed by i, where q i represents the charge number of the ion (e.g.q i = +1 for a monovalent cation like Na + ).
For conditions of dilute ionic congurations, this equation can be simplied by linearizing the exponential term in eqn (1) and leveraging the fact that the salt formula is electroneutral i:e: P i q i ¼ 0 . The resultant mathematical description predicts an exponentially decaying electrostatic potential, where z is the separation from the electrode interface, V ext is the applied potential, and l D f1= ffiffiffiffiffiffiffi ffi c ion p is the Debye screening length of the electrolyte solution.Gouy-Chapman (GC) theory, 21,22 as described in eqn (2), is applied widely when considering screening phenomena in the EDL, especially when they are of relevance to interpreting the results of electrochemical experiments.GC theory is particularly attractive because it advances the same simple physical picture as Debye-Hückel (DH) theory, 23 i.e., that a single tagged charge is surrounded by a diffuse neutralizing "cloud" of mobile ionic counter-charges.A single characteristic length scale emerges in both GC and DH, which is the so-called Debye length, l D .While physically intuitive, the theory for potential decay presented in eqn (2) involves a number of strong assumptions, many of which may be violated in electrochemically-relevant electrolyte systems. 7,8,24,25irst, GC theory assumes that the ions in the electrolyte can be modeled as point charges that occupy no volume.Second, GC theory assumes that the solvent environment can be accurately represented by a dielectric continuum with single dielectric permittivity, 3, neglecting both the nite size of solvent molecules, as well as possible correlations between their positions and orientations. 7,26While these assumptions may be accurate in some cases, they can result in qualitatively inaccurate predictions on the length scale that characterizes the EDL.For example, in 1 M monovalent aqueous electrolyte, the Debye length l D = 3 Å, which is comparable to the hydrated radius of a single solvated ion; at this scale, the intuitive physical picture of a diffuse ionic charge screening cloud of width l D becomes untenable. 24,26Additionally, GC theory does not model the strong molecular interactions between a planar electrode surface and solvent dipoles, which have been shown to exhibit strong orientational preferences within a few molecular layers of an electrode surface. 25,27,28ome of the issues associated with the second assumption are remedied by the so-called Stern correction, which posits that there is a molecular layer (the "Stern layer") of specically adsorbed ions at the electrode surface. 8,29Including the Stern correction in the theory developed in eqn ( 2) results in an electrostatic potential prole that decays linearly in the Stern layer, and then exponentially out to the bulk.Fig. 1C depicts schematic descriptions of potential decay as described by each of the theories discussed.

Dependence of the interfacial potential profile on ionic concentration
We study electrostatic screening in the EDL using molecular dynamics (MD) simulations of aqueous NaCl electrolytes at various different concentrations, conned between two Pt electrodes.Fig. 2A depicts a representative snapshot of the molecular simulation cell, with the z axis oriented normal to the planar electrode surface.In our simulation model, the electrodes are held at constant potential, as if they were connected to an external potentiostat.1][32] In this method, the partial charges on the electrode atoms are adjusted between each MD timestep in order to maintain a constant potential difference between the two electrodes.For the simulations presented in this section, the le electrode is held at a potential of V L = −V ext and the right electrode is held at V R = V ext , where V ext = 0.5 V, thus imposing an overall potential difference of DV = 1 V between the two electrodes.
The simulation box has dimensions 3.1 nm × 3.1 nm × 9.3 nm, and is periodically replicated in the directions lateral to the electrode surface.Previous studies have shown that the local dielectric constant of water approaches the bulk value within 35 Å of an interface. 33With a wall separation of 80 Å our simulation cell is large enough to host two non-overlapping EDLs at each of the electrodes, with a well-formed bulk in the central region.The number of water molecules and ions in the simulation cell vary with electrolyte concentration and are presented in Table 1 in the Methods section.The intermolecular potential, or "force eld" describing the water molecules is the standard TIP3P force eld, 34 whereas the ion and Pt atom interaction parameters are taken from studies reported in the literature. 35,36These parameters were chosen due to the excellent correlation of calculated and experimental ion hydration properties for this particular force eld. 36Ion solvation properties are known to play an important role in ion distribution preferences at interfaces. 10,37The results presented here depend on the details of the intermolecular interaction potentials, but the qualitative screening structures reported should be broadly conserved between force elds.
Our primary basis of comparison between continuum theory and atomistic simulation is the average electrostatic potential prole, i.e., the Poisson potential, 4. In this work, we represent this prole as an average over an equilibrium ensemble of single snapshot potentials.Specically, we dene the mean Poisson potential as, 4ðx; y; zÞ ¼ 1 N X N i¼1 4i ðx; y; zÞ; where 4i (x,y,z) is an instantaneous representation of the Poisson potential for the i th MD snapshot.Computing 4i (x,y,z) for a single conguration of a point-charge species requires rst dening a spatially continuous charge density prole, ri (x,y,z), on a regular lattice.We accomplish this with a proportional charge spreading scheme, which interpolates the charges on atomic centers to nearby points on the regular lattice.The details of this charge-spreading scheme are described in the Methods section.Aer dening ri , we compute 4i from the Poisson equation, V 2 4i = −r i /3 0 .The averages presented below are derived from equilibrated simulations of 10 ns with snapshots taken every 100 ps (i.e., N = 100).We dene the onedimensional potential prole as, where L x and L y are the length of the electrolyte region in the x and y directions, respectively.Dening the potential in this way enables straightforward comparison to the one-dimensional proles predicted from continuum theory and admits the analysis of microscopic potential uctuations via the statistics of 4(x,y,z).Fig. 2B depicts traces of the plane-averaged Poisson potential for several different electrolyte ion concentrations.The thin black lines represent the approximate planes of constant potential imposed in the MD simulation, where the voltages are pinned at V L = −0.5 V and V R = +0.5 V.It is apparent that the electrostatic screening structure is well established over the course of the simulations, with two distinct EDL regions conned within 10-15 Å of the electrode surfaces, and a bulk region with a at electrostatic potential prole, thus indicating the absence of static elds from the electrodes.Fig. 2C shows zoomed-in snapshots of the EDLs at the le and right electrodes, highlighting the electrostatic potential variation in these regions.Although the plane-averaged potential proles exhibit electrostatic potential decay from the electrode surface to the bulk, in line with continuum descriptions of electrostatic screening, they are strikingly dissimilar from the proles depicted in Fig. 1.First, the proles are markedly non-monotonic, exhibiting relatively large oscillations in the local electrostatic potential on the scale of V ext .Second, and perhaps most strikingly, the electrostatic potential proles show little variation over the entire range of ion concentrations studied here, instead of becoming more compact at higher concentrations, as predicted by GC theory.
The marked deviations from GC behavior in Fig. 2 raises the question: in what manner does GC-like screening manifest in the atomistic system, if at all?To answer this question, we attempt to separate the contributions of ions and water molecules to interfacial screening.We accomplish this by analyzing the screening prole of ions only, via the construction of the ionic screening function, which for the le EDL is dened as, where z L denotes the position of the le electrode, and r ion (z) is the xy-planeaveraged ionic charge density, obtained by restricting the proportional spreading procedure (described in Section 6) only to the atom-centered charges on the ionic species, and neglecting the charges of the solvent molecules entirely.
The screening function for the right electrode is dened analogously to that of eqn ( 5) but with an integral extending from z R to z R − z.Intuitively, the screening function S(z) tallies the amount of electrode charge that remains unscreened by the mobile ionic charges of a given distance into the bulk; it takes the value Q (elec) at the electrode surface (z = 0), and levels off upon reaching the bulk region of the simulation cell.We observe that S(z) z 0 for values of z beyond the EDL width.Fig. 3A depicts traces of the normalized ionic screening function S(z) h S L (z)/ S L (0) for various different ionic concentrations.According to GC screening theory, this normalized screening function should decay from unity to zero in an exponential manner, with associated length scale l D , the Debye length.The proles in Fig. 3A show signatures of the GC screening behavior; the proles are more diffuse at lower salt concentrations, and analysis of the length scale by exponential tting, as illustrated in Fig. 3C, shows that the decay length is roughly l D for each ion concentration.
Although the ionic statistics are in line with expectations from GC theory, the molecular dynamics simulations provide ample evidence that the solvent molecules play an important role in electrostatic screening.Fig. 3B shows traces of the bulk-normalized concentration of water molecules (solid lines) and sodium cations (dashed lines) in the simulation for the different ionic concentrations examined.At these concentrations, the rst density peak of the water molecules appears closer to the le electrode than the rst density peak of the cations.Additionally, the height of the peak, normalized to the bulk density of water, is roughly independent of the ion concentration, indicating that the water molecules are able to screen the electrode charge at distances closer than the typical cation approaches the electrode.

Deconstructing water's influence on the interfacial potential profile
In continuum theories, such as those based on the seminal work of Debye and Hückel, 23 the role of solvent in determining the electrostatic potential prole is reduced to that of a simple dielectric medium. 7,8Any notable variations in the shape of the potential are thus attributed to the spatial redistribution of mobile charge carriers.The results in the previous section reveal that water plays a more signicant role than ions do in shaping the potential prole within the EDL.To isolate this role, we consider simulations of neat water conned between two electrodes under varying applied potential.
The neat water simulations utilize a similar simulation size and setup to the system presented in Section 3, but use the SPC/E force eld for water and electrode atom force eld parameters for graphite rather than platinum, while still maintaining constant potential in the same way.The SPC/E force eld was used because it provides a dielectric constant closer to experiment than TIP3P.The switch to graphite was motivated by the observation that water chemisorbs onto platinum at positive potentials as demonstrated by ab initio simulations. 38,39ater chemisorption is expected to inuence the quantitative details of the interfacial potential prole, and this feature cannot be captured with classical force elds.However, water does not chemisorb onto graphene at the potentials used in our simulations. 40Therefore, we expect our simulations will at least qualitatively capture the reorientation of interfacial water, which is expected to contribute to the inner-layer capacitance of the interface. 41We carried out simulations with the electrodes held at constant potentials of V ext = 0.00, 0.25, 0.50, and 1.00 V, imposing an overall potential difference of DV = 0.00, 0.50, 1.00, and 2.00 V.Additional details are provided in the Methods section.Again, we analyze these simulations by computing the mean Poisson potential prole, 4(z).
The potential prole computed for unbiased conditions, i.e., V ext = 0, is plotted in Fig. 4A.This potential prole exhibits oscillations near the electrode, similar to those appearing in Fig. 2. The bulk potential for neat water at V ext = 0 is at with a value of 4 bulk z −0.55 V. We note that this interfacial potential drop and the oscillations are not accounted for in the treatment of water as a simple dielectric continuum.
A potential prole computed under a 2 V bias (V ext = 1 V) is plotted in Fig. 4B.This prole exhibits similar features to those of the unbiased prole with the addition of a nite slope in the bulk.This slope implies the presence of a static electric eld, thereby indicating incomplete screening of the applied electrode potential by water, as expected for a neutral solvent.This partial screening is reminiscent of the effect of a uniform dielectric on the eld between a parallelplate capacitor, as we discuss further in Section 4.3.However, the observed slope in the potential differs from the expectations of dielectric continuum theory.
Taken together, the potential proles plotted in Fig. 4 reveal three primary contributions that water makes in shaping the interfacial potential prole.These contributions are: (1) pronounced oscillations over molecular length scales, (2) a roughly 0.5 V drop in potential over the rst 1 nm of the interface at both electrodes, leading to a bulk-level potential that is not at the midpoint of the electrode potentials, (3) a reduction in the electric eld within the bulk under applied electrode potential.In the following subsections, we discuss the molecular origins of each of these contributions and how they relate to the dielectric properties of the water-electrode interface.

Molecular scale oscillations in the interfacial potential prole
Both ab initio and classical molecular dynamics simulations have revealed molecular scale oscillations in the electrostatic potential prole of solid-liquid interfaces. 6,27,42These oscillations are attributed to the consequences of molecular layering at the electrode surface, which serves to break translational symmetry.An increased density (relative to the bulk) of water molecules at the contact plane of the electrode exclude the adjacent plane, resulting in a subsequent decrease in density.This phenomenon, which has an analog in the oscillation of a radial distribution function, results in the emergence of well-dened hydration layers, as illustrated in Fig. 5.
The oscillating water density prole leads to a similarly oscillating charge density eld.When subject to the Poisson equation, this oscillating charge prole naturally results in an oscillating potential prole.The length scale of oscillation is determined by the molecular size and can be modeled with classical density functional theory. 43Oscillations like this, but persisting well beyond l D , emerge in ionic liquids. 12Since these oscillations have been previously well studied using MD simulation, we refrain from elaborating on them further herein.

Water's interfacial potential drop
One striking feature of the unbiased potential prole plotted in Fig. 4 is the potential drop between the neutral electrodes and the bulk liquid.This potential drop is also evident in the biased system plotted in Fig. 4B and in the aqueous electrolyte systems plotted in Fig. 2. Analysis of molecular dynamics simulation data reveals that this effect originates from anisotropy in the orientations of water molecules in the rst hydration layer.This anisotropy is intrinsic to the waterelectrode interactions and thus symmetric between the le and right electrodes.The anisotropy is apparent in the distribution of molecular orientations plotted in Fig. 5, which reveals a preference for interfacial water molecules to direct their hydrogens toward the electrode (away from the bulk).This orientation prevails because it provides favorable coulomb interactions between the image charges in the electrode and the partial positive charges of the hydrogen atoms.
The orientational bias of interfacial water molecules leads to a charge density wave that is net charge neutral, as illustrated in Fig. 5.According to the Poisson equation, a neutral density wave of this form (negative charge oriented toward Fig. 5 (A) and (B) The orientational distribution function for the angle, q, made between the water dipole vector and the electrode surface normal, as illustrated schematically above.The differently shaded lines correspond to populations of water molecules a given distance, Dz, from the electrode surface.The dashed grey line is the distribution corresponding to bulk water.(C) and (D) The charge density profile is computed for the population of water molecules at the interface of the electrode.Blue and red lines correspond to charge densities at V ext = 0.0 V and V ext = 1.0 V, respectively.Snapshots of the interface over the same horizontal axis scale are included to establish a sense of molecular lengths.
increasing z) yields a nite potential drop.At the opposite electrode, this charge density wave is mirrored leading to an equal and opposite potential rise.The quantitative details of this effect depend sensitively on the water-metal interactions as well as the distribution of charge within the water molecule.The magnitude of this effect is thus expected to depend sensitively on the choice of force eld and on the simulation conditions.
To understand this effect, we consider a simple model of dipolar solvent polarization at constant potential boundaries.This model includes a onedimensional charge density prole extending along the z coordinate, as illustrated in Fig. 6.The charge density prole is described on a lattice with lattice spacing ', representing the approximate radius of a water molecule.Under unbiased conditions with neat water, both the electrodes and bulk liquid have charge densities of r = 0. We model the charge density wave associated with solvent polarization with a discrete charge density wave, r(−L) = r w and r(−L + 1) = −r w , with a symmetric contribution at the other electrode, r(L) = r w and r(L − 1) = −r w .
The potential prole that results from the lattice charge density wave is plotted in Fig. 6B.The potential prole of this lattice charge density is sigmoidal, becoming constant at z ¼ 2' with a bulk potential value of 4 bulk ¼ Àr w ' 2 =3 0 .Taking ' ¼ 2:3 Å, the polarization density required to achieve a value of 4 bulk = −0.55V is r w = 0.57e nm −3 .
If the inherent orientation of water molecules at the interface is inverted, i.e., with hydrogen atoms pointed away from the electrode, then the resulting charge density wave results in a bulk potential that has a higher value than that of the neutral electrode.

Partial screening of external electric elds within the bulk
A dielectric medium, such as liquid water, has the general effect of reducing the strength of electric elds originating externally, such as from extended charged surfaces, or internally, such as from charge solutes.If a medium with dielectric constant 3 is exposed to an external electric eld of magnitude E 0 , then the eld within the bulk medium is E bulk = E 0 /3.Liquid water has a large dielectric constant, 3 w z 80, owing to its large molecular dipole and can thus signicantly reduce the strength of external elds.
However, when liquid water is located between parallel electrodes held at constant potential, the dielectric effect simultaneously amplies and reduces the elds originating from the electrode.The reduction arises due to the polarization of solvent dipoles and the amplication arises due to the concomitant charging of the electrodes in order to maintain their potential.For a uniform dielectric continuum, these effects exactly cancel, yielding no net reduction in the bulk electric eld.To maintain this, the electrode charge density q el = 3q 0 , where q 0 is the charge density required to maintain potential in vacuum (3 = 1).
The electric elds in our simulations of neat water under constant potential bias are lower than expectations based on a constant potential dielectric continuum, yet higher than expectations based on constant charge parallel plate capacitors.This is illustrated in Fig. 4. Fig. 7 contains a plot of E bulk vs. applied potential, indicating that the partial screening of the external electrode elds is a linear function of applied potential.
Analysis of our MD simulation data indicates the inuence of liquid water on the electrode charge density is similar to that of a standard dielectric medium, yet not identical to that of a homogeneous dielectric continuum.The presence of water results in an amplication of electrode charge density relative to the case where the electrodes are separated by a vacuum, and the strength of the electric eld from this amplication is signicantly reduced in the bulk relative to that of bare electrodes.However, the strength of each of these effects implies a different value of the dielectric constant.Furthermore, neither effect is consistent with the dielectric constant of this model of liquid water.
The partial screening of the electric elds from the electrodes can be quanti-ed in terms of an effective dielectric constant, 3 (E) eff , which is dened in relation to the slopes of the potential originating from the bare electrodes and the potential within bulk liquid.Specically, E bulk = E el /3 (E) eff , where E el = q el /3 0 is the eld originating from an innite plate with surface charge density q el and E bulk is the slope of 4(z) within the bulk, as indicated in Fig. 4. The value we compute for this effective dielectric constant is 3 (E) eff z 60, and roughly independent of applied potential.Specically, 3 (E) eff = 59.8 when DV ext = 0.5, 3 (E) eff = 61.5 when DV ext = 1.0, and 3 (E) eff = 60.1 when DV ext = 2.0.This effective dielectric constant is similar in magnitude to that reported in the literature for this water model, 3 SPC/E = 73. 44he amplication of electrode charge due to water can be quantied in terms of a different effective dielectric constant, 3 (q) eff , which is dened in relation to the charge density of the electrodes separated by liquid water and separated by vacuum.Specically, q el = 3 (q) eff q 0 , where q 0 is the charge density to maintain the potential in vacuum.The value we compute for this effective dielectric constant is 3 (q) eff z 24, and also roughly independent of applied potential.Specically, 3 (q) eff = 24.4 when DV ext = 0.5, 3 (q) eff = 24.3 when DV ext = 1.0, and 3 (q) eff = 23.4 when DV ext = 2.0.This effective dielectric is signicantly smaller than that of the liquid.
We assert that these seemingly incommensurate observations arise from deviations from dielectric continuum theory due to molecular effects at the water interface.As a dielectric continuum, water is assumed to be spatially uniform and everywhere charge neutral.Both of these assumptions break down at the charged electrode-water interface.Excluded volume effects limit the plane of closest approach for water molecules and orientational anisotropy due to the external electric eld leading to narrow planes of charge buildup with equal magnitude and opposite sign at either electrode.Fig. 8A illustrates the concept and Fig. 8B demonstrates that there is indeed a positive net charge buildup at the negative electrode.There is an analogous negative buildup at the positive electrode.The position of this charge plane is displaced from the electrode surface by a nite distance, ' ¼ 2:3 Å, and its magnitude scales with V ext .
This physical picture implies a simple one-dimensional model of the charge density eld of an electrochemical cell with two opposing polar solvent-electrode interfaces, as illustrated in Fig. 9.In this model, one electrode is located at position z = −L and held at potential V = −V ext and the other electrode is located at position z = L and held at potential V = +V ext .The charge density of the electrodes are represented by Gaussian distributions centered at z = −L and z = L: Àq el e ÀðzþLÞ 2 =2s 2 þ q el e ÀðzÀLÞ 2 =2s 2 ; and the water polarization layers are represented by similar distributions displaced into the bulk by a distance ', where s is atomic in scale.The full charge density eld is given by a sum of these two contributions, as illustrated in Fig. 9.The value of q el is determined by the constant potential condition that 4(L) − 4(−L) = 2V ext and the value of dq determines water's interfacial dielectric response.
The potential differences between the electrode and the center of the bulk can be computed directly from the Poisson equation.In this model, by construction, 4(0) − 4(−L) = 4(L) − 4(0) = V ext .In the limit that the Gaussian distributions are innitely narrow, this expression is simply, where L ′ = L/3 0 and ' 0 ¼ '=3 0 .The '/0 limit of this expression corresponds to the standard dielectric continuum picture.In this limit, q el − dq = V ext /L ′ = E 0 .Noting that E 0 = q 0 /3 0 = 3 −1 q el /3 0 , we see that dq = q el (1 − 3 −1 ).Ð z 0 r w ðzÞdz, illustrating that at the negative electrode, there is a positive net charge from water.There is an equal and opposite net positive charge at the adjacent positive electrode (not plotted).We denote the distance, ' between the first peaks of r el and r w which contributes to the effective dielectric response of the water-electrode interface.
The deviations from dielectric continuum theory that we highlighted above (i.e. that 3 (E) eff z 60 and 3 (q) eff z 24) arise from the case where 's0.In this case, eqn (9) can be manipulated to yield, thus indicating that the degree of interfacial solvent polarization is directly proportional to that of the electrode.The proportionality constant, a, depends on ', L, and the effective dielectric constant, 3 (q) eff .Eqn (10) provides a basis for understanding the unexpectedly low value of 3 (q) eff = 24.We note that dq is related to water's bulk dielectric constant via a simple parallel plate capacitor model, where 3 w is water's bulk dielectric constant.Substituting dq from this expression into eqn (9) yields, This expression reveals that the relationship between 3 (q) eff and 3 w is mediated by the ratio '=L, with 3 (q) eff = 3 w in the '/0 limit.In our system, where L = 8 nm, 3 w = 73, and 3 (q) eff = 24, we nd 'z2:3 Å, which is consistent with the results plotted in Fig. 8.The physical picture that this model advances is like that of the Stern layer.Some fraction of the applied potential drops due to water's interfacial dielectric response in the '-wide region between the electrode and the rst solvent plane.The magnitude of this drop is effectively equal to dV ¼ 4ðÀLÞ À 4ðÀðL À 'ÞÞzq e1 '=3 0 .The remaining potential drop, DV GC = V ext − dV, is thus what remains to be screened by the migration of ions (e.g., following Gouy-Chapman theory).In a pure water system, this quantity is, This expression has the feature that for a macroscopic system, i.e., ' ( L, DV GC = V ext .In other words, dV / 0 as '=L/0.

The role of ions in amplifying water's interfacial dielectric response
The implication of the above analysis is that the interfacial potential drop due to water (i.e., dV) is only appreciable in nanoscale systems.This conclusion is specic to a pure solvent system, where there are no mobile charge carriers to participate in screening.In this case, any unscreened potential drop (e.g., DV GC ), must extend across the entire length of the system.Potentials dropped over macroscopic length scales require only small electrode charges.In reality, neat water contains dilute concentrations of "water ion", i.e., H 3 O + and OH − , that can contribute to attenuating potential drops to microscopic scales (l D ∼ 1 mm in water at neutral pH).When this is the case, the analysis above must be adjusted to account for the inuence of ionic screening.Assuming the rst hydration layer excludes ions (so Gouy-Chapman-like screening occurs at the le electrode starting from z ¼ ÀL þ ') and that the potential is fully attenuated over the distanceð' þ lÞ; 4ðÀðL À ð' þ lÞÞÞ À 4ðÀLÞ ¼ 4ðLÞ À 4ðL À ð' þ lÞÞ ¼ V ext , and thus 4ðL À ð' þ lÞÞ À 4ðÀðL À ð' þ lÞÞÞ ¼ 0. Eqn (9) can thus be revised as, where l ′ = l/3 0 .It can be shown that with ionic screening q el = 3 0 V ext /l, whereas without ionic screening q el = 3 0 V ext /(3L).With this expression, the interfacial potential drop, dV ¼ V ext ð'=ð' þ lÞÞ and likewise, implying that the amount of potential that must be screened by ions is a decreasing function of ionic strength (increasing function of l).If ' ¼ 2:3 Å, then according to this expression, DV GC z 0.63V ext for a 0.6 M solution (l D = 3.96 Å) and DV GC z 0.52V ext for a 1.5 M solution (l D = 2.47 Å).

Conclusions
The predictions of Gouy-Chapman theory are based on the assumption that water is a homogeneous dielectric continuum.In this manuscript, we have used allatom molecular dynamics simulation to demonstrate that this assumption breaks down in multiple important ways.Water molecules take up physical space which leads to the appearance of molecular layering at solid-liquid interfaces and an associated displacement of the solvent polarization layer away from the electrode into the bulk.In addition, the molecular orientations of water molecules at the interface are anisotropic, owing to the specic details of non-spherically symmetric water-metal interactions.Together, these effects play a dominant role in shaping the interfacial potential prole, thus obscuring the comparatively minor effects of classical Debye-Hückel-like ionic screening.

Simulation details for systems with varying ionic concentration
Simulations at varying concentrations included the number of ions and water molecules indicated in Table 1.

Simulation details for neat water systems
The dimensions of the simulation cell were 2.8 nm × 2.9 nm × 9.3 nm.Water was modeled with the SPC/E water model 45 and the number of water molecules was 2042.The charges on the electrode atoms were allowed to uctuate to maintain constant potential with the ELECTRODE package 46 in LAMMPS. 47For each V ext = 0.00, 0.25, 0.50, and 1.00 V, we generated 100 trajectories each 1 ns in length.From each of these trajectories, the last 0.8 ns were used for data analysis.This provided an aggregate total of 80 ns of simulation data for each V ext .

Computational details for Poisson potential computation
Given a molecular simulation trajectory, we would like to devise a numerical scheme to determine the Poisson potential at any point in the simulation volume.
The Poisson potential 4 is dened by the Poisson equation, where r is the free charge density eld.The geometry under consideration has two periodic dimensions, denoted x and y, and one closed dimension z, which is the coordinate normal to the planar electrodes.The appropriate boundary conditions are, 4(0,y,z) = 4(L x ,y,z) 4(x,0,z) = 4(x,L y ,z) where V le and V right are the applied potentials on the le and right electrodes, situated at z = 0 and z = L z , respectively.To start developing a numerical scheme, we can discretize eqn ( 16) on a threedimensional rectangular grid.Given a specication of the number of grid points N = (N x ,N y ,N z ), we can dene a vector of grid spacings D h Note that the z-coordinate has a slightly different grid spacing because we would like to impose Dirichlet boundary conditions in this coordinate, requiring an extra boundary point.Now, grid points can be indexed by an index tuple n {0,.,N x − 1} × {0,.,N y − 1} × {0,.,N z }.The spatial location of a grid point with index tuple n is simply r n = D$n.Eqn ( 16) can be naturally discretized using a second-order Laplacian stencil.Under a row-major indexing scheme for the coordinates, this produces a sparse representation of the Laplacian operator.
Solving eqn ( 16), discretized on a grid, is relatively straightforward if we have a way to evaluate the free charge density eld, r, on the grid points.However, particles in a molecular simulation are, in general, not situated on a uniform grid.Hence, we need a scheme for interpolating a non-uniform charge density eld onto a uniform grid of points.Formally, for a particle with index k carrying charge q k localized at position r k = (x k ,y k ,z k ), we identify eight points bounding the voxel containing the particle.The index tuples of these eight points can be computed using the following equations, Now, we assign each point bounding the voxel a weight, The charge density on all grid points is computed by incrementing the charge density on each point i ˛1,.,8 bounding the voxel containing particle k by its charge weight q k w k i , and repeating for all particles in the simulation.

Poisson-Boltzmann theory for electrode charge
We can adapt Poisson-Boltzmann theory to generate a potential prole for the system architecture summarized in Fig. 9.To begin, consider the case where ' ¼ 0, i.e., traditional continuum approximation.In this case, linearized Poisson-Boltzmann theory yields potential prole of, According to this model, the potential due to the electrode is given by, which implies (through a parallel-plate capacitor model) a surface charge density for the le electrode of, q el = −V ext /l D = −(dq + q ion ), where q ion is the integrated excess ionic density prole in the interfacial screening layer.As the second equality suggests, charge neutrality necessitating that the charge contribution arising from the solution (i.e., dq + q ion ) must exactly counter the electrode charge.
We can adapt the equations above to describe the case where ' .0. In this case, the electrode charge is separated by the onset of the Poisson-Boltzmann screening layer (eqn (40)) by a distance '.The potential dropped over this distance is dV ¼ q el ', leaving the remaining potential DV GC to be dropped by solvent/ions.According to eqn (42), the charge in the solution will be, dq + q ion = −DV GC /l = −q el , (43)   where again, the second equality arises due to charge neutrality.Now, noting that DV GC ¼ V ext À dV ¼ V ext À q el ', we nd that, q el ¼ ðV ext À dV Þ=l ¼ ðV ext À q el 'Þ=l and thus q el ¼ V ext =ð' þ lÞ.With this,

Fig. 1
Fig. 1 Schematic electrical potential profiles in the EDL predicted by (A) Helmholtz-Perrin theory, (B) GC theory, and (C) GC theory with the Stern correction, where the dashed section represents the linear decay in the Stern layer.

Fig. 2
Fig. 2 Poisson potential computations from molecular dynamics simulations.(A) Representative snapshot of a molecular dynamics simulation cell, depicting an aqueous 1 M NaCl electrolyte between two Pt electrodes.The left electrode is the cathode (held at V L = −0.5 V), and the right electrode is the anode (held at V R = +0.5 V).By convention, the z axis runs perpendicular to the electrode surfaces.(B) Traces of the equilibrium-averaged, plane-averaged Poisson potential 4(z), estimated from molecular dynamics simulations run at various electrolyte salt concentrations.(C) Traces of 4(z) zoomed in on the cathode and anode, highlighting oscillatory and concentration independent short distance screening behavior.

Fig. 3
Fig. 3 Signatures of ionic and dipolar screening behavior.(A) Traces of the normalized screening function, S(z)/S(0) at various electrolyte concentrations, reflecting concentration dependent ionic screening behavior.(B) Traces of the local concentration of water molecules (solid lines) and Na + cations (dashed lines), normalized by their respective bulk concentrations, at various electrolyte concentrations.(C) Empirical screening length, l derived by fitting S(z)/ S(0) to a decay function exp(−(z − 8 Å)/l).The Gouy-Chapman prediction, where l = l D , is plotted as a solid black line.

Fig. 4
Fig.4The potential profiles for neat water with and without applied potential.(A) The Poisson potential, 4(z), for neat liquid water between neutral electrodes (V ext = 0).The calculated bulk potential, 4 bulk z −0.55 V. (B) Poisson potential for neat water with an applied electrode potential difference of 2 V (V ext = 1 V).The bulk electric field, E bulk , is the slope of the calculated potential (black line) evaluated with the bulk region.For reference, schematic potentials with slopes E 0 = V ext /L and E 0 /3 w are plotted in blue and green, respectively.

Fig. 6 A
Fig. 6 A schematic model of the influence of a symmetric charge wave on the electrostatic potential profile.(A) A piecewise representation of a water-like interfacial charge wave localized within 2' of the interface.(B) The Poisson potential arising from the density profile in (A) has a systematic drop in the bulk, analogous to that observed in the analysis of atomistic simulations.

Fig. 7 A
Fig. 7 A plot of the electric field in the bulk liquid at different values of the applied electrode potential DV = 2V ext .The plot also includes an indication of the predictions from dielectric continuum theory between surfaces of fixed potential, E 0 , and fixed charge, E 0 / 3 w .

Fig. 8 (
Fig.8(A) A schematic depiction of an interface between an electrode and a neutral dipolar liquid.In the top half and bottom half of the schematic, the net dipole orientation is rendered as arrows and partial charges, respectively.At a negative electrode surface, the orientation of the first hydration layer results in a narrow plane of net positive charge at the interface, subsequent planes of alternating charge cancel out as the degree of molecular layering diminishes.(B) Charge density profiles derived from simulation data of neat water at the left electrode with potential −1 V. Purple and green lines indicate the normalized charge density contributions of electrode atoms r el (z) and water molecules r w (z), respectively.The blue line represents the cumulative water charge density, r * w ðzÞ ¼

Fig. 9 A
Fig.9A model for understanding the interfacial dielectric response of a polar liquid at constant potential boundaries.(A) A depiction of the charge density field which combines the electrode surface charge, with amplitude q el and the water polarization surface charge, with amplitude dq.(B) The potential profile derived from the Poisson equation in the limit that the Gaussian distributions narrow to delta functions.The slope of the potential in the bulk region is that of a parallel plate capacitor with charge density r el − dr.

n k 1 =) n k 5 =) n k 6 =) n k 7 =) n k 8 =
(Qx k /L x S,Py k /L y R,Pz k /L z R) (23)n k 2 = (Px k /L x R,Py k /L y R,Pz k /L z R)(24)n k 3 = (Qx k /L x S,Qy k /L y S,Pz k /L z R)(25)n k 4 = (Px k /L x R,Qy k /L y S,Pz k /L z R) (26(Qx k /L x S,Py k /L y R,Qz k /L z S) (27(Px k /L x R,Py k /L y R,Qz k /L z S) (28(Qx k /L x S,Qy k /L y S,Qz k /L z S) (29(Px k /L x R,Qy k /L y S,Qz k /L z S),(30)where P$R and Q$S represent the integer oor and integer ceiling functions, respectively.Along each dimension d ˛(x, y, z), the particle position partitions the line segment connecting two adjacent grid points into two segments, one of length' that D ðdÞ ¼ ' ðdÞ Y þ 'ðdÞ[ , due to the properties of the ceiling and oor functions.For notational convenience, dene,

Table 1
Number of atoms contained in each simulation at the specified concentrations