Non-Reciprocal Interactions Induced by Water in Confinement

Water mediates electrostatic interactions via the orientation of its dipoles around ions, molecules, and interfaces. This induced water polarization consequently influences multiple phenomena. In particular, water polarization modulated by nanoconfinement affects ion adsorption and transport, biomolecular \hl{self-assembly}, and surface chemical reactions. Therefore, it is of paramount importance to understand how water-mediated interactions change at the \hl{nanoscale. Here we show that near the graphene surface anion-cation interactions do not obey the translational and isotropic symmetries of Coulomb's law.} We identify a new property, referred to as non-reciprocity, which describes the non-equivalent and directional interaction between two oppositely charged ions near the confining surface when their positions with respect to the interface are exchanged. Specifically, upon exchange of the two ions' positions along the surface normal direction the interaction energy changes by about 5$k_BT$. In both cases, confinement enhances the attraction between two oppositely charged ions near the graphene surface, while intercalation of one ion into the graphene layers shifts the interaction to repulsive. While the water permittivity in confinement is different from that in bulk, the effects observed here via molecular dynamics simulations and X-ray reflectivity experiments cannot be accounted for by current permittivity models. Our work shows that the water structure is not enough to infer electrostatic interactions near interfaces.


I. INTRODUCTION
Water and ions, two basic components of living organisms and ubiquitous in numerous natural and technological processes, are strongly interconnected [1][2][3][4]. Water mediates interactions between ions, molecules, macromolecules, and interfaces through the polarization of its electric dipole moment. At the same time, charged species induce polarization of other molecules, surfaces, and the medium itself [5]. This leads to complex interactions, the details of which are not yet fully understood. In particular, the capability of a medium to be polarized and the resulting attenuation of electrostatic forces are typically quantified by the dielectric permittivity. This dielectric permittivity normalized by that of the vacuum (ε 0 ), called the relative permittivity (ε r ), is frequently assumed to be a distinctive property of a material and is often termed as the dielectric constant. Such a denomination, however, is not applicable in nanomaterials, and especially for water in confinement and near interfaces where the properties differ from the bulk [6][7][8][9].
The underlying motivation to quantify a permittivity is to understand molecular interactions. Hence, considerable efforts have been made to determine the in- * e-mail: fenter@anl.gov † e-mail: m-olvera@northwestern.edu terfacial permittivity [10][11][12][13][14] of confined water. These studies have shown that the interfacial permittivity is anisotropic [12,13] and is significantly diminished within 1 nm from an interface [11]. However, the connection between this interfacial dielectric permittivity and the molecular interactions is not established. Here we investigate ion-ion effective interactions near a solid-water interface and the relationship between these interactions and the interfacial water polarization. We show that the symmetries of ion-ion interactions in bulk break down near an interface. Consequently, the description of molecular interactions via a dielectric function is not viable near an interface.
We select graphene to study ion-ion interactions near an interface. Graphene is technologically relevant because it possess extraordinary properties such as partial wetting transparency, electrical conductivity, and mechanical strength [15,16]. The graphene-water interface is of interest for water desalination [17,18], for electrochemical energy storage [19] and harvesting [20], as a trans-electrode membrane to characterize biomolecules [21], and in many other applications. Extensive studies of ions at air-water [22,23] and at graphene-waterair interfaces [24] have discussed ions size, polarizability, and solvation energy on their adsorption to hydrophobic surfaces. Here, we demonstrate that water polarization plays a central role on ion specific interactions and adsorption.
First, we investigate the water polarization near a arXiv:2009.00496v1 [cond-mat.mtrl-sci] 1 Sep 2020 graphene surface in the absence of any free charges and observe the changes of the interfacial water polarization in the presence of a nearby ion. Second, we compare the molecular interfacial water structure predicted by the simulation with that determined by X-ray reflectivity measurements. Third, we evaluate the effective ion-ion interactions near the graphene surface both in the plane and along the surface normal and assess the effect of increasing the confinement. We explain our results in terms of the interfacial water polarization rather than the dielectric constant. We demonstrate that the symmetry breaking of ion-ion interactions is a general phenomena which is present in symmetric non-polarizable ion models and in polarizable models of water, ions, and grephene. Finally, we analyze our results in the framework of the prevailing continuum theories of electrostatics at interfaces and we demonstrate that continuum theories fail to capture a fundamental breaking of symmetries of ion-ion effective interactions near interfaces.

A. Molecular Simulation Models
To understand the water-mediated and confinementaltered interactions, we investigate ion-graphene and ionion interactions for a pair of non-polarizable test ions (monovalent anion and cation of equal size) placed between two uncharged non-polarizable graphene surfaces (see Figure 1a). These assumptions allow us to show that the phenomena uncovered here are not attributed to specific properties of ions, such as size and polarizability. The generality of our findings is further confirmed by considering polarizable models of ions, water, and graphene (see Appendix A). The ions are investigated in different configurations with respect to the graphene surface. The simulation box consists of two graphene surfaces separated by a water layer of thickness L w ; each surface is formed by four graphene layers, i = 0, . . . , 3, with an inter-layer separation of 0.358 nm and 1008 carbon atoms per layer. The simulation box size is 5.065 nm × 5.104 nm × L z , where the box length along the z-direction, L z , is adjusted according to the number of water molecules, N w ; N w = 8060 and 2233 for L w = 9.8 and 3 nm (L z = 12.28 and 5.63 nm), respectively. The water is simulated using the SPC/E model [25], and the graphene parameters are taken from the OPLS-AA force field [26]. The ions' Lennard-Jones parameters are σ = 0.333 nm and = 1.16 × 10 −2 kJ/mol for both ions. Unless explicitly stated, the simulations are performed using the 3dc-method with slab correction to mimic 2d periodic boundary conditions in the x and y directions [27]. The system temperature is maintained at T = 298 K using a Nosé-Hoover thermostat.
To enable a one-to-one comparison between X-ray reflectivity experiments and MD simulations, a separate model was built with no added ions in the system and The response of the water polarization to electrical fields is investigated by placing two oppositely charged ions in different configurations. Here the cation (blue sphere) is at about 0.5 nm from the left graphene surface while the anion (green sphere) is at z ≈ 8 nm; anions and cations have the same Lennard-Jones parameters (σ = 0.333 nm and = 1.16×10 −2 kJ/mol) and differ only in their valence. Lz is the simulation box length in the z direction, and hi with i = 1, 2, 3 are three regions of 0.5 nm employed to quantify the water polarization next to the graphene surface. (b) Irregular surface employed for a direct comparison between X-ray reflectivity experiments and MD simulations; the surface is made of four partial graphene layers (G0 -G3) with the relative surface area of each graphene layer determined by the experimental XR results. (c) Calculation scheme employed to calculate the polarization field around an ion nearby the water-graphene interface.
with a graphene surface composed of four layers with different partial coverages (Figure 1b). This mimics the experimental surface wherein the graphene sample is made of multiple incomplete layers. The coverage of each layer was set to the value determined by the experimental Xray reflectivity best fit structure as follows: G 0 contained 1176 carbon atoms with x × y dimensions 5.157 nm × 5.955 nm (1 monolayer, ML); G 1 contained 924 carbon atoms, x × y = 5.157 nm × 4.679 nm (0.78 ML); G 2 contained 336 carbon atoms, x × y = 5.157 nm × 1.701 nm (0.29 ML); and G 3 contained 84 carbon atoms, x × y = 5.157 nm × 0.425 nm (0.07 ML). These surface coverages are similar to those determined experimentally.
which includes the dipole moment per unit volume P 1 , the quadrupole moment P 2 , octupole moment P 3 , and all higher-order moments. In water, the dipole moment is the main contribution to the polarization. Hence, it can be expressed as where µ i is the dipole moment of the i-th water molecule at r i , and . . . represents the ensemble average. The SPC/E water dipole moment is µ 0 ≡|µ i |=0.0489 e·nm where e is the positive elementary charge. The calculation scheme in molecular dynamics (MD) simulations is depicted in Figure 1b. Operationally, the instantaneous polarization is calculated as p(r, z) = ∆m ∆V , where ∆m = i µ i includes the dipole moment (µ i ) of all the water molecules in the volume ∆V at (r, z). The mean polarization is calculated as P(r, z) = p(r, z) , performing the ensemble average over at least 10 4 independent configurations.

B. Potential of Mean Force
In an N -body system the mean force exerted on the ith particle is derived from the instantaneous forces from all the particles and is given by where W (n) (r n ) is the n-particle (n ≤ N ) potential of the mean force (F i ), ∇ i is with respect to the coordinates of the i-th particle (i ≤ n); V (r N ) = N i<j u(r ij ) is the systems potential energy and r N represents the 3N particles' coordinates; u(r ij ) is the particles' pair interaction energy; . . . (n+1) represents the ensemble average over n + 1 . . . N particles. The potential of mean force is related to the n-particle probability distribution function G (n) (r n ) by In our study n = 2 and particles 1 and 2 represent the two ions. The potential of mean force depends on the separation distance between the two graphene surfaces, L w , and the ions' positions with respect to the graphene surfaces, r 2 and r 1 ; d = r 2 − r 1 is the ions' relative position and d = |r 2 − r 1 |. We investigate the following cases: 1. W (z) = W (r 1 , r 2 ) is the interaction of one ion, say 1, and the graphene surface along the surface normal direction, where |r 1 | = z 1 1.5 nm, |r 2 | → L w ; L w 10 nm. The interaction of ion 2 with the surface is calculated by exchanging the ions' position.
in the plane parallel to the graphene surface when z 1 = z 2 ≈ 0.28 nm, is the ion-ion interaction in bulk (i.e., with no the graphene surfaces).
To calculate the potential of mean force (PMF) we use the umbrella sampling method [29] which is based on Eq. (4). The technique consists in performing biased sampling by fixing the ion of interest at designated positions along the reaction coordinate ξ = z, d. The test ion is fixed using a harmonic potential given as u umbrella (ξ) = −k(ξ −ξ i ) 2 where {ξ, i = 1, . . . , M } is a set of M equilibrium positions and k = 2000 kJ/(mol nm 2 ) is the spring constant. The separation distance between two contiguousξ i andξ i+1 equilibrium position is 0.03 nm, approximately. The system is simulated over at least 20 ns to generate a distribution functions around for eachξ i . The biased potential of mean force is related to the particle's distribution functions within the windows by Eq. (4). The unbiased PMF, the biased PMF, the biased probability distribution function, and the external potential are related in an exact way [30]. The WHAM algorithm is employed to construct the unbiased PMF [31].

C. Continuum Theory of Electrostatics
In the classical theory of electrostatics, an interface is modeled by taking into account two media of dielectric constants ε 1 and ε 2 , respectively [28]. The dielectric mismatch between the two media leads to the following boundary conditions for an electric field E passing through the interface where E ⊥ and E || represent the electric field components perpendicular and parallel to the interface, respectively, passing through medium 1 or 2 as indicated by the superscripts. The boundary conditions in Eq. (5) imply that the work needed to bring a particle of charge q from infinity to a distance z from an interface of planar geometry is given by where the dielectric mismatch is quantified by α = (ε 1 − ε 2 )/2(ε 1 + ε 2 ) [32,33].

D. Molecular theory of permittivity
In the molecular theory of dielectrics, the interfacial region is described by means of a local permittivity ε [12] which is a function of the system's molecular parameters. The components of the permittivity tensor are calculated from the unperturbed interfacial water structure, i.e., in the absence of free charges within the interfacial region (See appendix B). In a system with slab geometry (see Figure 1a), the permittivity tensor is a diagonal matrix with components ε || (z) ≡ ε xx (z) = ε yy (z), and ε ⊥ (z) ≡ ε zz .
and for the inverse permittivity ε −1 ⊥ (z) in the direction perpendicular to the surface where || and ⊥ indicate, respectively, the parallel and perpendicular components of the total dipole moment M and the instantaneous polarization p(z), . . . means the ensemble average, β = 1/(k B T ), T is the absolute temperature, k B is the Boltzmann constant, ε 0 is the vacuum permittivity, and V is the water volume; (7) and (8) are applicable in systems where 3D periodic boundary conditions are assumed (x, y, and z directions). Ballenegger and Hansen [12] derived the expression for 2D periodic boundary (x and y directions) where D ⊥ = 1. The results from simulations using 3D and 2D boundary conditions are quantitatively different.

E. Experimental Methods
High-resolution X-ray reflectivity (XR) experiments provide a sensitive probe of the molecular structure of well-defined interfaces via a direct measurement of the electron density distribution [34] and have been employed to study the interfacial water structure at numerous planar surfaces [35][36][37][38]. Such measurements then serve as an indirect probe of the molecular dipole orientation and polarization at the interface. Specular XR measurements of the graphene/water interface, i.e., along the graphene (0001) surface normal direction, were carried out using a 3 mm × 10 mm sample of an epitaxial graphene thin film grown on a semi-insulating 6H-SiC(0001) substrate (EG/SiC) immersed in ultra-pure de-ionized water (DIW, nominal pH = 7). EG was grown by thermal decomposition of SiC under a partial pressure of Ar at 1500 The specular X-ray reflectivity (XR) is measured as a function of the perpendicular momentum transfer Q = |K f -Ki| where the incident Ki and reflected K f X-ray beams vary by the scattering angle 2θ.
• C according to methods previously described [39,40]. This approach produces incomplete graphene layers as a growth artifact. Nevertheless, the resulting sample with large-area graphene surfaces is better suited for XR measurements than mechanically exfoliated graphene, which produces micron sized flakes. Measurements were performed at beamline 33-ID-D of the Advanced Photon Source at Argonne National Laboratory. A photon energy of 14 keV (λ = 0.89Å) was used, and the X-ray beam cross section measured approximately ∼ 50 µm × 1 mm (h × w). The specular XR data probe the timeand laterally-averaged (in the xy-plane) interface structure as a function of the vertical momentum transfer, Q = 4π sin(2θ/2)/λ where 2θ is the X-ray scattering angle (see Figure 2). The sample area probed by the X-ray beam varies with the scattering angle and ranges from ∼ 1.58 mm 2 to ∼ 0.14 mm 2 for the present measurements.
The XR best fit structure was derived following a nonlinear least squares optimization procedure. In general, the electron density distribution of each atomic layer, j, along the substrate surface normal direction is modeled as a Gaussian with parameters to describe its position, z j , coverage, θ j , and thermal widths, u j . The total electron density profile is given by where Z j is the atomic number of the j th layer. The reflectivity signal R(Q) is related to the modulus squared of the Fourier transform of the electron density distribution and can be calculated for any set of model input parameters as Here, the proportionality indicates that experimental considerations such as the angle-dependent transmission of X-rays through the sample cell and surface roughness must be accounted for as extrinsic factors in the model (see Appendix C); f 0,j is the atomic scattering factor of each element Z j [41] and accounts for the Q-dependent decay of the X-ray scattering intensity resulting from the spatial distribution of electrons around an atomic core. We evaluate the accuracy of the model via the goodnessof-fit metric χ 2 where N is the number of data points, R exp (Q) is the experimentally measured reflectivity, and σ(Q) is the experimental uncertainty at Q. The model parameters are refined until χ 2 converges (χ 2 = 1 for a perfect fit within experimental uncertainties). The model for the present system consists of a semiinfinite bulk 6H-SiC(0001) substrate [42], one unit cell of SiC (i.e. six alternating C-Si layers) to describe surface, up to eight C layers with the density of 2D graphene, and a semi-infinite layered water model described by Magnussen et al. [43]. The interfacial SiC, graphene, and water structures are optimized according to Eqs. (10)- (12) while the bulk SiC structure is fixed.
Given a single adsorption surface (i.e., for a uniform, complete graphene layer), the layered water model includes a series of m Gaussians (m = 0, 1, 2... with the zeroth layer being closest to the adsorption surface) along the surface normal. The position z m and width u m of each layer are given by where z 0 is the height of the zeroth Gaussian relative to the substrate surface (i.e., the interfacial water height), d w is the distance between adjacent Gaussian peaks, u 0 is the width of the zeroth layer, andū is the width broadening of subsequent layers such that the density asymptotically approaches that of bulk water (ρ w = 330 e − /nm 3 , where e − is the negative elementary charge). We assume that the areal density of water in each Gaussian of this layered water distribution has the same coverage due to a lack of confinement in the lateral directions in our system. The coverage of each layer is then given by where A U C is the unit cell area of the SiC substrate and V w = 0.0299 nm 3 is the effective volume of a water molecule in bulk assuming spherical symmetry. As noted previously, the EG/SiC growth methodology used in this work leads to partial layers of graphene. Therefore, we modified the layered water model [38] to incorporate multiple graphene surfaces, G n (n = 0, 1, 2, ...), and assumed that water interacts in the same way with each graphene layer. Namely, above each exposed graphene surface exists the same intrinsic water structure according to Eqs. (14) and (15) but with a modulation due to the position z n , width u n , and coverage Θ n of the adsorbing graphene layer. All together, the layered water structure is described by where the graphene monolayer coverage on SiC is Θ ML = 3.147 carbon atoms per SiC unit cell. Because the least-squares fitting finds a local minimum in the parameter space, multiple structures consistent with the same XR data are possible. Therefore, we constrained the model parameters of the SiC and graphene surfaces based on previous chemically-resolved measurements of the air/EG/SiC interface [37]. Further details of the XR analysis and all best-fit parameter values are reported in Appendix C.

A. Interfacial water polarization
The polarization field of water near the graphene surface and in the absence of nearby test charges is shown in Figure 3a. Even if the surface is uncharged and there are no nearby ions and no external electric field, the water is polarized perpendicularly to the graphene surface. This intrinsic interfacial polarization is due to the water dipole moment alignment caused by the solid surface and is also observed when polarizable models of water and graphene are used in the simulations (see Appendix A). This preferential orientation decays in an oscillatory way to a random polarization in the bulk (z > 1 nm; see Fig. 3a). As a reference for ions in bulk we show the polarization field around a cation and around an anion far away (z > 2 nm) from the interface in Figures 3b and 3c, respectively. In this case, the test charges are sufficiently far from the graphene surface and do not perturb the interfacial properties. At this separation distance from the graphene surface the water polarization around the ions is spherically symmetric and decays without oscillations. However, the polarization field is not spherically symmetric when one ion is placed within ∼1 nm of the graphene surface (see Figures 3d and 3e for the nearby cation and anion, respectively) while the oppositely charged ion is kept away (z ≈ 8 nm). Moreover, the polarization response is asymmetric and unequal around the anion and cation. This result contrasts with the continuum theory of electrostatic where the polarization fields for positive and negative test ions are assumed to have the same magnitude and opposite direction.
The unequal responses are quantified by the potential of mean force (PMF), W (z), and is different for the anion and the cation (even if both have the same radius).   around the cation, and (c) around the anion; the arrows indicate the polarization orientation and the magnitude is given by the scale bar at the right; r is the distance from the ion center along the direction parallel the graphene surface (e || ) while z is the distance from the graphene surface in the perpendicular direction (e ⊥ ). The water polarization near the graphene surface is changed by the presence of a nearby (d) cation or (e) anion. The bottom panels show the ion potential of mean force profile to bring an ion from the bulk to the graphene surface from our molecular dynamics simulations Ws (purple line) and calculated via Eq. (6) form the continuum theory of electrostatics Wc (light-blue line).
For example, the cation's PMF at z ≈ 0.5 nm is approximately 0.5 k B T while the anion's PMF is approximately 0 k B T . Additionally, for z 0.5 nm the PMF becomes steeper for the anion than for the cation. The interfacial ion specificity has been attributed to the ionic polarizability, size, and valence. However, our results show that solely a change in sign gives rise to a pronounced ion specificity via the asymmetric water polarization response to the sign of the ions. In continuum electrostatic theory the change in the dielectric permittivity at the interface is modeled by the image charge method, which from water to graphene would assign an equally repulsive force for the negative and positive ions [28] (see light-blue line in Figs. 3d,e). This symmetry is also obtained in continuum models for electrolytes in confine-ment [33,44,45].

B. Interfacial Water Electron Density Distribution
Here, we focus on the details of the water structure in the absence of ions and the relationship between XR experiments and the simulation results. Analysis of the XR data ( Figure 4a) reveals a graphene/water electron density profile with four partial surfaces to which water adsorbs (Figure 4b). These include three layers of graphene, G 1 , G 2 , and G 3 with fractional layer coverages of 0.84 ML, 0.31 ML, and 0.1 ML, respectively, and a reconstructed carbon buffer layer, G 0 with complete coverage and which separates true 2D graphene (  Figure 1a). The water orientation is given by cos Ω, where Ω is the angle between a water molecule dipole moment (µ) and a unitary vector normal to the graphene surface (e ⊥ ); µ0 ≡|µ|=0.0489 e·nm is the dipole moment of the SPC/E water model, where e is the positive elementary charge; from left to right the three insets show representative configurations of water for cos(Ω) = −1, 0, and cos(Ω) 0, respectively. (e)-(g) The histograms of the instantaneous water polarization components H(pi) (i = x, y, z). Px and Py are always equal to zero. The polarization histograms in the z-direction are persistently narrower than in the x-and y-directions up to a water layer thickness of Lw ≈ 120 nm.
from the SiC substrate below (see the complete interface structure in Appendix C and Table I for details).
In general, the XR measurements and MD simulations reveal qualitatively equivalent water distributions adsorbed on the partial graphene surfaces (Figure 4b).
Both the XR-derived structure and the MD prediction show a weakly modulated water profile with density peaks that correlate with the locations of the graphene layers. The XR best-fit structure shows water adsorbed at ≈ 0.31 ± 0.03 nm above each exposed graphene surface while MD predicts the adsorbed water height to be 0.33 nm, both of which are consistent with a slightly hydrophobic interface [46,47]. The XR data analysis also reveals that the water may adsorb closer to the G 0 "buffer" layer than to the subsequent free-standing graphene layers by 0.01 ± 0.05 nm. Although this difference is smaller than the sensitivity of the XR measurement, it suggests that G 0 may be slightly less hydrophobic than graphene. This property of G 0 is supported by previous XR and ab initio MD (AIMD) results of water adsorbed on EG/SiC by Zhou et al. [38].
Zhou et al. found, via XR measurements, a significantly reduced water height of 0.23 nm on the G 0 buffer layer compared to a height of 0.32 nm above free-standing graphene, which suggested that the G 0 layer exhibits a hydrophilic character. In contrast, their AIMD results showed a ∼0.02 nm decrease in the water height above G 0 (z ∼0.31 nm) compared to that above freestanding graphene layers (z = 0.33 nm), which is in agreement with our XR results. Zhou et al. attributed the discrepancy in G 0 water height between their XR and AIMD results to surface defects in their EG/SiC sample, which likely also explains the observations in the present study. We discuss sources of the G 0 water height discrepancy in greater detail in Appendix C. Finally, MD predicts a broadening of the first hydration layer for thicker graphene regions (i.e., on G 3 where the total thickness of the graphene slab is ∼1 nm ), a phenomenon not observed in the XR results, with the root mean square (RMS) width of water on G 3 being most consistent with the experimental results (see Appendix C). Potential explanations of the subtle differences between our MD prediction and XR results may reflect different interactions between water and graphene of various thicknesses and finite size effects of the simulation compared to the relatively large area of the XR measurement (see Appendix C).
The intrinsic water structure is the density profile corresponding to a uniform graphene surface (see Figure 1a) extracted from XR measurements. Figure 4c shows a comparison of the intrinsic water density profiles from XR experiments and MD simulations. Both the XR and MD oxygen density profiles reveal a first hydration layer (z ≈0.3 nm) with a peak density that is more than twice that of the bulk. The density oscillations decay rapidly with a small secondary hydration layer at z ≈0.6 nm and a nearly bulk-like third hydration layer at z ≈1 nm in both XR and MD profiles. The atomic positions of hydrogen are calculated via MD simulations only due to the weak scattering of X-rays from hydrogen atoms. As such, the excellent agreement between experimental and computational oxygen distributions supports the accuracy of the proton distributions and the non-zero interfacial water polarization predicted by our MD simulations. The position of the oxygen layers relative to the graphene surface (validated by the experiment) are used to build the histograms of the interfacial water polarization which is employed to explain the electrostatic interfacial interactions (see below).
The location of the hydrogen density peak at the first hydration layer coincides with the oxygen density peak but is broader, suggesting that a fraction of the water dipole moments are oriented perpendicularly away from the graphene surface. This is partially driven by transient hydrogen bonds between water molecules of the first and second hydration layers (see left inset in Figure 4c) and the inability of water molecules to form hydrogen bonds with the graphene surface. The water orientation is quantified by the distribution function of its dipole orientation cos Ω with respect to the graphene surface normal e ⊥ (Figure 4d). The distribution reveals a first hydration layer in the region h 1 (z < 0.5 nm) with the majority of water dipole moments oriented parallel to the graphene surface (cos Ω = 0). However, the asymmetry in the distribution exhibits a slight preference of the water dipole moment to orient away from the graphene surface. The water molecules in the h 1 region form a hydrogen bond network which is a characteristic behavior of water near hydrophobic molecules [48] (right inset in Figure 4c). In the region defined at 0.5 nm< z ≤ 1 nm (h 2 ) the preferential orientation of water molecules is diminished (Figure 4c), and from z > 1nm all the water dipole moment orientations are equally probable (1 nm< z ≤ 1.5 nm).
We now examine in further detail the interfacial water polarization. In recent experiments it is investigated the out-of-plane permittivity (ε ⊥ ) of water confined between two flat surfaces of graphite and hexagonal boron nitride. The experiments reveal that ε ⊥ ∼ 2 within an interfacial layer of two or three molecules thick (referred to as the electrically dead layer) which is considered to have a "vanishingly small polarization" [11]. Here in contrast, we observe a non-zero persistent water polarization in the z-direction within the interfacial region. Figures 4e,f,g show the histograms of the instantaneous polarization H(p) within the regions h 1 , h 2 , and h 3 , respectively, defined to capture the different hydration layers. The average polarization in the z-direction is P z ≡ p z ≈ 0.03e/nm 2 in the h 1 -region (Figure 4e), P z ≈-0.004e/nm 2 in the h 2 -region (Figure 4f), and P z =0 in the h 3 -region ( Figure 4g) and for z 1 nm, while P x and P y are always equal to zero. In particular we highlight the non-vanishing polarization P z within the h 1 region which is fundamental to understand the electrostatic interactions near an interface (see next section).
The dipole fluctuations are known to be suppressed along the surface normal direction [49]. Here, the sup- pression is seen as a much narrower histogram for the polarization in the z-direction than in the x-and ydirections. The polarization suppression implies that every fluctuation of µ z is nearly canceled by an antiparallel component. The suppression is observed across the entire water region and up to a water layer thickness of L w ≈ 120 nm. Additionally, we investigated the liquidvapor interface by removing the graphene surfaces, and we found that the dipole moment fluctuations in the zdirection remain suppressed. This observation suggests that the suppression is a general result for the propagation of electrostatic fields across interfaces and is not related to the chemical structure of the confining surfaces.

C. Ion-Ion Interactions in Confinement
We look into a direct force analysis between an anion and a cation in bulk water, near the graphene surface, and between an ion in the water and an ion intercalated between the graphene layers. Ion intercalation is the common mechanism for Faradaic energy storage (e.g., lithium ion batteries with graphitic anodes) [50], while ion adsorption occurs in capacitive energy storage systems. Figures 5a-c show that confinement enhances the ion-ion interaction along the surface normal direction, which is more pronounced at L w = 9.8 nm than in the bulk and even more pronounced by decreasing the water layer thickness to L w =3 nm; the minimum of the PMF decreases by about 5k B T by decreasing L w from 9.8 nm to 3 nm. The enhancement of the ion-ion interaction near the graphene surface is consistent with observations near a hydrophobic surface [4]. Furthermore, the effective ion-ion interaction under confinement is non-reciprocal, i.e., the potential of mean force is not equivalent upon exchange of the ions' positions (See Figure 5a, b) where it is seen that the minimum of the PMF decreases by about 5k B T by exchanging the ions' positions. This nonreciprocal behavior of the anion-cation interactions is also observed when polarizable models of Na + Cl − ions, polarizable water, and polarizable graphene are used in the simulations (see Appendix A). The ion-ion interaction in the surface plane direction W || , however, is not significantly different from the ion-ion interaction in bulk (Figure 5c) likely because the ions are less confined in the xy-plane. Interestingly, the interaction along the surface normal direction changes from predominantly attractive when both ions are in the aqueous phase (Figure 5a,b) to completely repulsive when one ion is intercalated within the graphene interlayer space (Figure 5d). In an isotropic medium the ion-ion effective force depends only on the ions' separation distance, and the forces are equivalent upon exchange of the ions' positions. However, figures 5a-d reveal that in confinement the effective ion-ion interaction is directional and position dependent with respect to the graphene surface. The force non-reciprocity cannot be explained by current models of anisotropic dielectric permittivity (ε || and ε ⊥ ), which predict the same dielectric profile regardless of the ion valence and, hence, a symmetric force. We will discuss this point in detail in subsection III E.
The non-reciprocal behavior of the ion-ion interfacial interactions and the ion-sign specificity are due to the unsymmetrical response of the interfacial water polarization. The classical theory assumes surface polarization is the same for positive and negative charges whereas our results show the polarization field strongly depends on the sign of the charge. The potential of mean force W c from the MD simulations takes into account the ion-ion direct interaction and the interaction through the surrounding water molecules, which includes the work to reorient the interfacial polarization. The interfacial polarization when the cation is placed on the graphene sur-face and when the anion is on the surface (exchanged configuration) are shown in Figures 6a and 6b, respectively. These polarization fields give rise to different ionion interactions. The polarization around two ions close to each other is canceled at a short distance ( 1.5 nm) which is consistent with the observation around the polar groups of proteins [51]. When the two charges are separated at longer distances, however, the polarization propagates over the whole separation distance ( Figure  6c). The polarization in confinement is enhanced even at larger confinement distances. Figures 6d,e show the difference in the water polarization at the middle distance between the two ions separated by about 9.5 nm in confinement and in bulk, respectively. The magnitude of polarization in confinement is approximately twice and shows less deviation from the perpendicular orientation than the polarization in bulk.

D. Results from continuum electrostatics theories
We calculate the ion-ion Coulombic interaction energy as a function of the distance between the two ions arranged in the three configurations in Figures 5a-c using the image charge method (see subsection II C). We assume that the water is a continuum background with the permittivity of ε 1 = ε w = 70 and the graphene is a continuum medium of permittivity of ε 2 = ε g = 2. The image charge method cannot predict the oscillations found in the potential of mean force in Figure 5 when surface polarization effects are included. Instead, a purely repulsive interaction is observed assuming a dielectric mismatch at the interface of α = 0.97. Figure 7 shows the ion-ion interaction energy (from the image charge method) when the cation is placed inside the graphene surface for different combinations of the dielectric constants of the two media forming the interface. It is interesting to see that the continuum theory captures the change from attractive to repulsive for sufficiently large dielectric mismatches. The repulsive behavior is in agreement with the results from umbrella sampling simulations when one ion is intercalated between the graphene layers (see Figure 5d), however, the interaction energy obtained from the image charge method is the same when the ions' positions are exchanged.

E. Results from the molecular theory of permittivity
We investigated the interfacial dielectric permittivity employing the molecular theory described in subsection II D. The expressions of the permittivity in the direction parallel to the graphene surface (ε || (z)) and the inverse permittivity in the perpendicular direction (ε −1 ⊥ (z)) are given by Eqs. (7) and (8), respectively. The ultimate goal of quantifying the interfacial permittivity is to determine the electrostatic interaction between charged species near  the interface. We estimate the interfacial ion-ion interaction by means of the Coulomb's potential employing a position dependent permittivity ε(r) given by Eq. (8). We see that the interaction potential W C significantly overestimates the interaction with respect to the potential of mean force from umbrella sampling simulations W ⊥ (see Figure 8). As we mentioned above, the poten-tial of mean force from umbrella sampling simulations W ⊥ is different by exchanging the positions of the two particles (non-reciprocal) whereas here we see that the ion-ion interaction calculated from the permittivity profile from Eq. (8) is the same independently of the ions arrangement (see Figure 8). This shows that the interaction between two charged particles can not be described solely in terms of a permittivity function.

IV. CONCLUSIONS
Our results show that in confinement the ion-ion effective interaction is directional and position dependent with respect to the graphene surface due to a non-zero persistent interfacial water polarization. We identify this property as non-reciprocity, which describes the directional dependent interactions and non-equivalent change of interactions upon exchange of the ions' positions near a confining surface. The non-reciprocity implies that ionion interactions at the interface do not obey the isotropic and translational symmetries of Coulomb's law and are observed in both polarizable and non-polarizable models. This phenomenon contrasts with the ion-ion interactions in an isotropic medium (bulk) where the force depends only on the ions' size and separation distance, is not directional, and is equivalent by exchange of the ions' positions. Traditionally, ion specificity is attributed to the internal ion polarization including polarizability associated with ion size. Here, we find that the water polarization plays a central role in the behavior of ions near interfaces. Namely, the water polarization around ions and near the interface alters electrostatic interactions, leading to non-equivalent ion-interface interactions upon exchange of the ion charge sign even if the ions have equal size. This non-symmetrical water polarization affects the understanding of ion-differentiation mechanisms such as ion selectivity and ion specificity. The agreement between XR experimental measurements and MD simulations of polarizable and non-polarizable models supports the presence of the interfacial water polarization. Current models based on the anisotropic dielectric permittivity do not account for non-reciprocal ion-ion interactions and confinement effects. Molecular interactions near interfaces and in confinement are related to a variety of processes including chemical reactions, adsorption, and transport phenomena, and these processes in turn are relevant in many technological and scientific problems. The insights gained here need to be considered in the understanding of processes based on asymmetric ionic adsorption and interactions at interfaces. The induced polarization of atoms results from the deformation of the electronic cloud due to a local electric field and it is known to be relevant in some interfacial phenomena [52]. At first approximation, the main contribution of the induced polarization is from the induced dipole moment (µ ind ) of atoms, ions and molecules due to the local electric field E'(r), given by where α is the polarizability, which is a distinctive property for each atom and ion. To account for polarization effects, we employ the SWM4-NDP polarizable water model [53] and the corresponding force field parameters for the polarizable Na + and Cl − ions [54]. Polarizable graphene is modeled by including the polarizability derived from DFT calculations [55] in our graphene model. The water molecule has a fixed HOH geometry bearing two positive charges at the hydrogen centers and the negative charge is placed at a fixed distance from the oxygen center, shifted along the molecule's axis of symmetry. The induced polarization is based on the Drude particle model of charge q D , attached to the center of the polarizable atoms (or ions) through a harmonic potential. The Drude particle's charge is balanced by the positive charge of the core such that, q C + q D = 0 for neutral atoms and q C + q D = q ion for ions, where q ion is the ionic charge. The spring constant k D , the polarizability α, and the charge of the Drude particle q D , are related by The polarizabilities employed in our calculations are α g = 1.139Å 3 for graphene carbon atoms, α w = 0.978 A 3 for water oxygen atoms, α Na + = 0.157Å 3 for Na + and, α Cl − = 3.969Å 3 Cl − . Integration of the equations of motion is performed by means of the extended Lagrangian algorithm [56] which consists in assigning a small mass to the Drude particles m D taken from the atomic masses, in such a way that the mass of the core is m i − m D . The Drude particles are simulated at a much smaller temperature than the whole system to fulfill the Born-Oppenheimer minimum energy condition. The Drude particles mass is m D = 0.4 g/mol. A dual Nosé-Hoover thermostat is used to maintain the Drude particles temperature at T D = 1K and the system at T = 298 K. We simulated a set of systems similar to that described in subsection II A, namely, a water layer of thickness L w ≈ 9.8 nm, formed by N w = 8060 water molecules, and confined between two graphene surfaces formed by four graphene layers each (see Figure 1a). A time step of 1 fs is employed for integration of the equations of motion of the systems without ions. When ions are present the time step is 0.2 fs. The simulation protocols are similar to those described in Section II. Consideration of polarizability of the water molecules implies an additional dipole moment contribution to the permanent dipole moment.   First we examine the water structure near the graphene surface considering four different combinations, namely, a) polarizable water and polarizable graphene, b) polarizable water and non-polarizable graphene, c) non- FIG. 11. Interaction between polarizable Na+ and Cl − ions at the water-graphene interface. The potential of mean force (W ⊥ ) as a function of the ion-ion separation distance d along the graphene surface normal direction. The two ions (Pin and Pout) are placed at different heights z above the graphene surface at the same x and y coordinates (see Figure 5a,b). The ion positions are exchanged such that in (a) Pin =Na+ and Pout =Cl − while in (b) Pin =Cl − and Pout =Na + . The ions' polarizabilities are α Na + = 0.157Å 3 and α Cl − = 3.969 A 3 for Na + and Cl − , respectively. The water polarizability is αw = 0.978Å 3 and the graphene polarizability is αg = 1.139Å 3 . The water layer is confined between two graphene surfaces separated by Lw = 9.8 nm.
polarizable water and polarizable graphene, and d) a) non-polarizable water and non-polarizable graphene. The density profiles as a function of the perpendicular distance to the graphene surface are shown in Figure  9. The four combinations predict qualitatively similar profiles exhibiting peaks at the same location, approximately. Our results are consistent with previous studies which show that the graphene polarizability does not significantly affect the water density profile [55,57]. Here, the polarizable water model predicts higher peaks the non-polarizable water. Analysis of the polarization field shows a persistent interfacial polarization qualitatively similar in the four combinations (see Figure 10).
In Figure 11 we analyze the effective interaction between two ions (Na + and Cl − ) in water nearby the graphene surface. All the components in the system (water, ions, and graphene) are polarizable. The ions are placed at different heights above the graphene surface aligned perpendicularly to the graphene surface in a similar way as in Figures 5a,b. The potential of mean force (PMF) is calculated in two cases: 1) the Na + ion is placed on the graphene surface while the Cl − ion is at different heights from the graphene surface and 2) the ions' positions are exchanged. The PMF profiles are qualitatively different from those for-non polarizable ions (see Figure 5a,b for L w = 9.8 nm). For example, the profiles exhibit less pronounced oscillations when Na + is on the graphene surface (purple line in Figure 11) than for the non-polarizable ions discussed Figures 5a,b. The difference in energy when the ion positions are exchanged, at the minimum of the PMF, is 5 k B T in the fully polarizable system (see Figure 11). This is similar in both magnitude and sign to that observed for the nonpolarizable system (see Figure 5a,b for L w = 9.8 nm). Hence, the non-reciprocal behavior of the ion-ion interactions is present in both systems with polarizable and non-polarizable atoms and ions. An electric field E applied in a dielectric material induces a polarization due to a separation of the bound charges in the material (atomic nuclei and their electrons). The macroscopic field D is called the electric displacement and is given as where P is the polarization density. Linear response theory assumes that for a weak applied electric field the, induced electric polarization is proportional to the magnitude of the applied field where χ is the susceptibility tensor, which is related to the dielectric permittivity tensor ε(r) by I is the identity matrix. D is then expressed as D = ε 0 ε · E (B4) Using index notation in Eq. (B4) in Cartesian coordinates can be expressed as In simple materials the dielectric behavior is isotropic, leading to a diagonal dielectric tensor with three equal components. Under nanoscale confinement the dielectric response is anisotropic, and the components of the dielectric tensor are not equal. In the next part we present the main steps to derive the expression for the dielectric permittivity for a non-uniform medium in one dimension. Our treatment follows closely the work by Feller and Stern [13], Ballenegger and Hansen [12], and Roland Netz' group [14]. The basic idea is to compute the response of a dielectric medium to a static, external, and uniform electric field E, combining descriptions from statistical mechanics and continuum electrostatics. By combining both approaches, we are able to derive the expression of the local permittivity.
Let ∆E(r) be the change in the mean local electric field inside the dielectric. When the external electric field E is turned on, the change is due to both the external field itself and the dipoles within the medium. E 0 (r) and P 0 (r) are, respectively, the mean electric field and the mean local polarization with no applied external field (E 0 (r) is zero if P 0 (r)). In the linear response regime the change in the local polarization ∆P(r) and the change in the total electric field ∆E(r) are related by where χ(r) is the local susceptibility tensor which is related to the local dielectric permittivity tensor ε(r) as The expression for ∆P(r) is derived from statistical mechanics while the expression for ∆E(r) is from macroscopic electrostatics.

Microscopic description
We consider a classical system in a microstate Γ described by the Hamiltonian H(Γ). In general, the dipole moment m changes from point to point within the dielectric. Hence, the instantaneous polarization density at r, p(r), (also known as electric polarization, or simply polarization) is given by The system's total dipole M is given by The mean polarization is expressed as Let H and H be, respectively, the system's Hamiltonian when there is no applied external electric field and when the electric field is turned on. In linear response theory H can be expressed as H = H − M · E, and the change in the mean polarization as In the weak electric field regime we linearize the above expression to get where we have switched to express the components α = x, y, z of the vector ∆P(r), and the statistical average is performed in the zero electric field regime. Eq. B12 involves the correlation between a fluctuation in the local polarization density m(r) and a fluctuation in the global dipole moment M of the system. To calculate the permittivity tensor via Eq. (B6) it is necessary to derive an expression for ∆E. The calculation of ∆E is performed by considering the dipolar contributions from every molecule and from each image cell in a system where periodic boundary conditions are assumed in the x, y, and z directions. By doing so, Feller and Stern arrived at the following expressions [13] for the dielectric permittivity profile in the parallel (ε || (z)) and perpendicular (ε ⊥ (z)) directions to the graphene surface given by Eqs. (7) and (8), respectively.
Appendix C: Experimental details

X-ray Reflectivity Background and Model
The specular XR signal R(Q) is directly related to the laterally-averaged real-space electron density profile ρ(z) via Fourier Transform as is the angle-dependent transmission of Xrays through the sample cell, B(Q) accounts for the effects of surface roughness [58], r e = 2.82 × 10 −5Å is the classical electron radius, and A U C is the unit cell area of the substrate (SiC in the present case); Q is the vertical momentum transfer as in Figure 2a. Due to the loss of phase information inherent in the XR measurement, it is not possible to determine ρ(z) by inverse FT. Instead, it is determined by optimizing a model wherein the j th atomic layer along the z-direction is represented by a Gaussian with fitting parameters of position, width, and coverage as described in the main text Eq. (10). The parameter values and their uncertainties are determined via non-linear least squares fitting following the Levenberg-Marquardt algorithm until the level of agreement (Eq. (12)) between the calculated model-dependent reflectivity, R(Q) calc , and the experimental data, R(Q) exp , converges. Figure 12 shows an AFM height scan of the EG/SiC sample in air. It reveals a smooth sample surface (average surface roughness = 38 pm) with lateral terrace widths of ∼ 1µm. We also see several layers of graphene on each terrace with a mean thickness of about 6Å and a maximum thickness of about 12Å (Figure 12, top). Due to the limited resolution of the AFM, it's chemical insensitivity, and the fact that it does not see the layers below the surface, we cannot determine the number of graphene layers precisely from the AFM. Nevertheless, this qualitative picture provides an additional reference point against which to evaluate the XR best fit structure (Figures 4a,b and 12, bottom). The optimized XR model parameters (z, u, and Θ) for each atomic layer are summarized in Table I, and the complete best fit interface structure of de-ionized water (DIW) on EG/SiC is shown in Figure 12. We partially constrained the surface SiC and G 0 parameters in the XR analysis based on the work of Emery et al. [37]. That study combined the chemical and structural sensitivity of X-ray photoelectron spectroscopy (XPS), X-ray standing waves (XSW), and XR to rigorously characterize the structure of the air/EG/SiC interface and address longstanding debates about the buffer layer, G 0 , between SiC and epitaxial graphene. The XPS/XSW measurements were consistent with a carbon-rich buffer layer composed of two chemically-distinct and partially overlapping layers, S1 with sp 2 hybridization and S2 with sp 3 hybridization and bonded to Si atoms of the substrate [59][60][61][62][63]. The results ruled out a proposed Si adatom model [36,64]. Although Emery et al. found that the EG/SiC interface was largely identical across all samples studied (UHV- sharp spikes in the height profile indicate the step edges. Bottom: resolution-broadened best fit electron density profile from the CTR data with χ 2 = 1.6 ( Figure 2a in main text), including SiC substrate results and the partial water layers adsorbed on each exposed surface; u broad = (u 2 i + u 2 res ) 1/2 where ures = 0.55/Qmax. grown vs. Ar-grown with 1.3 or 1.7 ML graphene), we allowed for the structure parameters in this work to vary up to 10% from their results. In general, our results are in agreement with those previously reported.

a. G0 Buffer Layer
We identified a G 0 layer with a mean height of 2.31 ± 0.02Å above the SiC surface and FWHM of approximately 0.72Å. Similarly, Emery et al. found a mean G 0 height of ∼2.3Å with a G 0 FWHM of ∼0.8-0.9Å. Others have reported G 0 heights of ∼ 2.5Å [60,64]. We note that although our XR measurement could not sufficiently resolve the two G 0 subpeaks S1 and S2, we were unable to obtain a good fit to the data using a single G 0 peak, likely due to the asymmetry of this layer. The resolution of the XR measurement is given by r = π/Q max ; Q max = 4.911 A −1 resulting in r ≈ 0.64Å for the current measurement, which is larger than the S1-to-S2 separation we identified of 0.42 ± 0.05Å but smaller than the FWHM of the combined G 0 layer. The S1-to-S2 peak separation is in agreement with the 0.35Å spacing reported by Emery et al. We found that the spacing between S1 and S2 was conserved throughout the fitting iterations as both layers moved together with respect to the SiC surface, lending support for the shape of the buffer layer. Moreover, the distance from S2 to the topmost Si layer of the substrate was 1.97 ± 0.05Å, in agreement with earlier reports [37,60]. We find a G 0 coverage equivalent to 1.18 ± 0.03 ML, in contrast to previous results that identified a layer with essentially the density of graphene [37,38,63]. Attempts were made to constrain the G 0 density to that of a single graphene layer, but such a density was always found to be inconsistent with the data. The excess carbon density we identified in the G 0 layer may account for a surface oxide species not included in our model [37,[65][66][67][68]. Emery et al. identified via XPS the presence of SiO x in several EG/SiC samples, which they were also unable to accurately model in the XR data analysis. They estimated the oxygen coverage to vary from 2 O/nm 2 in an Ar-grown sample to 6 O/nm 2 in a UHV-grown sample. The excess carbon density of the G 0 layer in our Ar-grown sample can equivalently be attributed to an oxygen content of ∼5 O/nm 2 .

b. SiC Surface
We identify a partially depleted SiC surface, consistent with the thermal desorption of Si during graphene growth [69][70][71]. We find that Si was depleted down to the fourth surface layer of the SiC while the C layers within SiC were not depleted. The topmost Si layer was displaced away from the bulk and toward the G 0 buffer layer. The coverage of the topmost Si layer is consistent with the coverage of the S2 layer of G 0 (though we report a large uncertainty on the S2 coverage, see Table I). We were unable to rigorously quantify the S2 layer coverage due to the large covariance of this parameter with the S1 coverage as a result of the limited resolution of the XR measurement discussed previously. Therefore, we fixed this parameter after several fitting iterations where it was converging to values consistent with the amount of surface Si depletion. That is, the coverages of the two layers indicate a one-to-one bonding between dangling Si atoms and sp 3 -hybridized carbons, consistent with the proposed growth mechanisms and previous reports [37,[69][70][71].

c. Epitaxial Graphene
We find a graphene film structure that is consistent with the AFM in Figure 12. Namely, we identify three 2D graphene layers, G 1 -G 3 , above the G 0 buffer layer for a total of 1.25 ± 0.07 ML of epitaxial graphene. The mean graphene layer spacing is 3.40 ± 0.07Å, in agreement with the known value [72], and G 1 is located 3.5 ± 0.04Å above G 0 , in agreement with previous reports [37,62,63].

d. Adsorbed Water
The overall water structure results are discussed in the main text, and the best fit parameter values are given in Table I. Here, we focus on similarities and differences between our results and a those from a previous XR study of the adsorbed water structure on graphene [38] and discuss the implications for the intrinsic interactions of water with graphene.
We find a water height above the free standing graphene layers (G 1 -G 3 ) of ≈ 3.1Å, in agreement with the results of Zhou et al. [38]. However, we identify significant differences in the water structure above the G 0 buffer layer. We find that water adsorbs closer to the buffer layer at ∼ 3Å above G 0 , but the height uncertainty overlaps with the water height above free standing graphene and indicates that the buffer layer is also weakly hydrophobic. We note that although we assumed a single water model above all "graphene" surfaces in our analysis (including G 0 ), the different height for the buffer layer results from calculating the distance between water and the weighted average position of the S1 and S2 peaks, a condition that was not imposed during the XR data analysis. In contrast, Zhou et al. found that water adsorbs at a height of ≈ 2.33Å above the buffer layer, indicating a hydrophilic character. They report water contact angle (WCA) measurements that support a more hydrophilic G 0 and which are equivalent to the WCA on bare SiC [73]. The WCA measurements showed a linearly increasing trend with graphene layer thickness, which may initially suggest that the buffer layer and subsequent graphene layers display wetting transparency on SiC. However, the extent of graphene's wetting transparency is still contested [16,[73][74][75], and it is unclear if the G 0 layer with its sp 2 and sp 3 character would exhibit similar wetting transparency properties. Instead, the evidence presented by Zhou et al. suggest that the G 0 water adsorption is more in line with a high concentration of defects on their samples.
Zhou et al. used UHV-grown EG/SiC samples, which are known to possess a greater amount of defects than graphene grown in a furnace in an Ar atmosphere (as was done for the sample studied in this work), as shown by their AFM images and other reports [76]. Based on the results of Emery et al. wherein UHV-grown and Ar-grown EG/SiC were found to have equivalent interface structures [37], which are consistent with our own EG/SiC interface structure, we would not expect the growth methodology to substantially contribute to the differences observed between our and Zhou's G 0 -water distance. However, the quality of the sample depends on the vacuum level and any pre-treatments of the SiC to remove oxides [65,66,68,77]. No pre-treatments were reported by Zhou et al. They also reported Raman data with significant D and D+D' peaks, which result from edge and other defect states [73,[78][79][80][81]. In fact, it has been shown that the introduction of such defect peaks upon oxygen plasma etching of EG/SiC is associated with a decrease in WCA [73]. Finally, Zhou et al. report MD and ab initio MD (AIMD) simulations of defect-free surfaces that predict water heights above G 0 consistent with that observed above free-standing graphene, and in agreement with our MD results. The AIMD simulations included effects of the SiC substrate and the corrugation of the buffer layer, but found only a ∼0.2Å decrease in the adsorbed water height above G 0 compared to freestanding graphene. Only upon inclusion of Si vacancies and -OH defects were they able to simulate a G 0 -water height of 2.33Å. We conclude that while we can reasonably expect water to adsorb more closely to the buffer layer than to the subsequent graphene layers as a result of the SiC substrate and corrugated surface, we expect the effect to be minor in the absence of substantial defects. For a defect-free surface, the water structures above the buffer layer and subsequent graphene layers are very similar and can be described well by a single water model given certain resolution limits of the XR measurement. Figure 13 shows adsorbed water profiles isolated from each exposed graphene layer (see Figure 2d in the main text for the total density profile). Although the XR and MD results generally agree, a noticeable broadening of the first hydration layer is predicted by MD above each partial graphene surface. The width is consistent with the XR result by the third layer (G 2 ). A graphene slab of four uniform and complete layers produces the same water structure as is shown in Figure 13 for the layer G 0 . The origins of these discrepancies should be discussed in the context of the limitations of both the MD and XR approaches:  (1) The G 0 layer of the MD simulation is modeled as a 2D graphene layer due to the absence of a SiC substrate. However, the partial sp 2 and sp 3 character of G 0 on SiC alters the band structure [82] and may affect the buffer layer's interaction with water. Indeed, this is consistent with our XR results and previous AIMD predictions [38], though the effect is weak. Previous MD simulations of unsupported graphene sheets suggest via water contact angle calculations that the graphene hydrophobicity decreases with increasing layer thickness [75] and becomes graphite-like by the third layer [83]. This agrees with our MD simulations wherein the adsorbed water peak height is conserved, but the broadening reveals a small increase in the number of water molecules that are able to adsorb closer to the graphene sheet.

Agreement and Differences with Simulation
(2)The XR-derived result is limited by the complexity of the XR analysis model, as described in detail above. We used a single intrinsic water model for each graphene surface despite evidence that on SiC thicker graphene regions are more hydrophobic than thinner regions [84]. We note that this is the opposite of MD predictions, suggesting that while the effect of the SiC substrate and corrugated buffer layer are relatively weak in magnitude, they may significantly affect the chemistry of the surface. Accounting for different hydrophobicity would introduce additional parameters to an already complex model and is unnecessary given that the changes are small. Using a single water model essentially gives an optimized structure that captures the average behavior of water adsorption on an imperfect graphene surface. In addition, the non-linear least-squares fitting of the data can produce multiple equally viable structures (i.e., with equivalently good χ 2 ). This point emphasizes the importance of including evidence from other experimental methodologies such as XPS/XSW [37] to refine and constrain models. Even with these caveats, the qualitative agreement between the XR and MD structures, paired with the high-confidence in the XR result (given a χ 2 of 1.6 where a perfect fit would have a χ 2 of 1), provides a consistent picture of water adsorption on graphene.