Origin of dielectric polarization suppression in confined water from first principles

It has long been known that the dielectric constant of confined water should be different from that in bulk. Recent experiments have shown that it is vanishingly small, however the origin of the phenomenon remains unclear. Here we used ab initio molecular dynamics simulations (AIMD) and AIMD-trained machine-learning potentials to understand water's structure and electronic properties underpinning this effect. For the graphene and hexagonal boron-nitride substrates considered, we find that it originates in the spontaneous anti-parallel alignment of the water dipoles in the first two water layers near the solid interface. The interfacial layers exhibit net ferroelectric ordering, resulting in an overall anti-ferroelectric arrangement of confined water. Together with constrained hydrogen-bonding orientations, this leads to much reduced out-of-plane polarization. Furthermore, we directly contrast AIMD and simple classical force-field simulations, revealing important differences. This work offers insight into a property of water that is critical in modulating surface forces, the electric-double-layer formation and molecular solvation, and shows a way to compute it.


S1 Simulation details
All relevant information about the simulations performed in this study are given in Table S1 below.
Table S1: Details of the AIMD FFMD and NNP simulations.For each slit thickness h, defined as the separation between the two confining materials, we provide the cell dimensions L x/y/z , the number of wall atoms N wall and water molecules N water , the equilibration period t eq and the total simulation length t sim .Slit

S2 Dielectric constant calculations S2.A Polarization computation
For all AIMD estimates of the dielectric constant reported in the main text, we computed the polarization from the charge density of our AIMD simulations.This assumes that the polarization for z = 0 is equal to 0, given that all atoms are far away.To validate this approach, we computed the total dipole moment of the simulation box for various configurations over a 500 fs simulation at an applied electric displacement of 1 Vnm −1 and h = 1.935 nm [S].The dipole moments were computed every 20 fs. Figure S1 contrasts the values obtained from the charge density as described in the main text, the Berry phase and the Wannier centers.It shows that P (z = 0) = 0 is a very good assumption and allows the computational cost of the Wannier centers to be avoided.The error made is of the order of 3.10 −3 atomic units for the dipole moment.Deriving the error on the polarization and using Eq.( 4) of the main text, this leads to an error on the inverse of the dielectric constant on the order of 10 −3 , thus validating our choice for computing the polarization.

S2.B Subtracting the surface contribution
In the AIMD simulations, the confining surfaces of the slits (graphene and hBN) contribute to the global polarization response due to the extension of their electronic clouds that partially overlaps with the one of water.In experiments, the confining material contributes also to the dielectric response, and its contribution is removed by performing a differential measurement [1].Thus, we used a similar approach here, where we corrected the calculated values by subtracting the local response of the system in the absence of water, that is composed only of the confining surfaces (see snapshots in Fig. S1a-c).First, we extracted the total local polarization response ∆P tot z from the total system.We then extracted the polarization response of the confining surfaces ∆P CS z from a single calculation without water.This enabled us to obtain an approximation of the water response as ∆P H2O z = ∆P tot z − ∆P CS z .By combining it with the expression of the local dielectric constant, we get The result of this procedure is reported in Fig. S2 below, showing the dielectric constant profile obtained calculating a single configuration for the case of graphene slits (Fig. S2, profiles on the left) and hBN slits (Fig. S2, profiles on the right).It shows that the contribution of the confining material to the dielectric response has been replaced by the one of vacuum.

S2.C Data convergence
To check the equilibration of the simulations, we analyzed the time evolution of the dipole moment of the simulation boxes along with its running average, shown in Fig. S3a below for water inside a graphene slit for the M system.The convergence of the total dipole moment was found to be very fast, as reported in a previous study on bulk water [2]: equilibration and convergence is faster at finite D than at finite E. We calculated the dielectric constant of the water layers reported in Fig. 1b in the main text from the running average of the dielectric constant, shown in Fig. S3b below, by taking the average of ϵ −1 ⊥ in the long-time limit, corresponding to the last 20 ps of each simulation, except for the S [0.911 nm] system where the calculations were made on the last 10 ps.For the L [3.00 nm] system, we used the last 60 ps, as the response of the second layer requires more statistics to converge.Even if Fig. S3b shows that the dielectric constant is not fully converged for the 1.935 nm [M] slit size, we have used it to compute the dielectric constant as it is at the limit of what can currently be reached with AIMD simulations.The error was estimated from the fluctuations of the running value of ϵ −1 ⊥ around this average, and propagated to obtain the error in ϵ ⊥ .
We summarized the obtained dielectric constants of the first and second water layers (L1 and L2, respectively) for graphene and hBN slits in Table S2 below.For graphene, we also computed the dielectric constant ϵ C,P CD ⊥ by assuming that the charge distribution corresponds to the one of the SPCE model, i.e. a point charge distribution Table S2: Dielectric constant of interfacial water layers.Calculated dielectric constants of the interfacial water layers, L1 and L2, reported in Fig. 1b in the main text, for both graphene (first column) and hBN slits (second column) and various water slab thickness using AIMD simulations.The third column computes the dielectric constant from the AIMD configurations by assuming that the charge distribution corresponds to the one of the SPCE model, i.e. a point charge distribution.The error is obtained from the fluctuations of the running average in the long-time limit.(PCD), in order to estimate the relative proportion of orientation and electronic effects.We found that the obtained values of ϵ ⊥ for graphene and hBN slits were equal within the error bar, except for the monolayer configuration which exhibits a slightly lower value for graphene.We obtained larger uncertainties in the dielectric constant for the second layer, L2, which we attribute to two factors.First, it arises from the higher error sensitivity of these calculations for higher values of ϵ ⊥ , that is, for a given uncertainty in ϵ −1 ⊥ , the uncertainty in ϵ ⊥ increases with increasing ϵ ⊥ .Secondly, this may arise from the more dynamical character of the system with larger slit thickness, with more molecules exchanging between the layers that slow down the convergence.The comparison with the dielectric constant obtained from the point charge distributions shows that most of the difference between the FFMD and AIMD calculations comes from the difference in the molecular orientation rather than electronic effects.
To understand the impact of the choice of DFT functional, we repeated the DFT simulations for the 1.935 nm [M] graphene slit using the BLYP functionnal with the D3 correction.These calculations yielded a dielectric constant for the first interfacial water layer (L1) of 2.096 ± 0.067, a value in agreement with and even lower than what we found in this study.We note however that the response of the L2 layer did not converge even with 120 ps simulations.

S2.D Impact of the dividing surface
The dividing surface, that is, the domain boundary over which the dielectric profile is integrated, is not well defined on the atomic scale due to the non-continuous nature of the interface between the confining layer and the water phase.We investigated the effect of the choice of the dividing surface by considering four positions, indicated in Fig. S4a (dashed lines): the confining surface (CS), the one used in the main text corresponding to the atomic position of the confining layer (e.g.carbon atoms for graphene); the electronic density crossover (EC), corresponding to the local density maximum between the carbon wall and the water phase; the full density (FD) surface, corresponding to the symmetric position of the EC surface with respect to the confining surface; and the jellium edge (JE) [3], corresponding to the extension of the electronic cloud of the confining material in a mean-field approach.
Figure S4b shows the calculated inverse dielectric constant for each position of the dividing surface against the experimental data as a function of the water slab thickness.The calculated values are sensitive to the dividing surface position, and this is due to the change in the domain volume, as the dipole moment remains largely unchanged.In particular, shifting the dividing surface from the atomic position towards the vacuum, as in the case of the FD surface, increases the volume at constant dipole moment, thus decreasing the resulting polarization, which in turn leads to a decrease in the obtained ϵ ⊥ .The opposite occurs when the dividing surface is shifted towards the liquid phase, as in the case of EC and JE surfaces: the polarization is increased with decreasing the volume, and in turn ϵ ⊥ increases up to infinity and even becomes negative.In this study we show that our calculations yield values of ϵ ⊥ in reasonably good agreement with the experimental values by simply taking the dividing surface at its physical position -the atomic position of the confining layer (CS) -and subtracting the electronic cloud contribution beyond the surface, as explained in section S2.A above, without any need of shifting the dividing surface to a different position.We note that shifting the dividing surface to the FD surface would lead to an even better agreement, but it would not reflect the experimental setup.Indeed, it would imply to compute the dielectric constant of the first water layer (L1) over a spatial domain larger than the experimental system, thus artificially reducing the dielectric constant.Importantly, the position of the dividing surface does not affect the calculation of the dielectric constant of the second water layer (L2), or the dipole moments, and the water molecular arrangement presented here.We also note that for all the dividing surfaces considered here, our AIMD calculations yielded smaller ϵ ⊥ than the values obtained from force-field simulations.

S3 Additional force-field simulations
We compared the out-of-plane dielectric constant, ϵ ⊥ , of confined water calculated from our AIMD simulations with values obtained from FFMD simulations with the aim of evaluating the performance of our simulations and the impact of different force-fields.Figure S5 shows different FFMD simulations, either reported in the literature or carried out in this study, of water confined in graphene slits as a function of the water slab thickness against the experimental data of Ref. [1].The water models used in the FFMD simulations reported here from the literature are SPCE [4] (used in [5,6,7]) and KKY [8] (used in [9]), a more complex model that allows intramolecular vibrations.Graphite and graphene are modeled by rigid carbon walls using the carbon-water two body interactions.The forcefield models are parametrized from various sources, namely, DFT calculations [9] using a Buckingham potential based on [10], properties of water droplets on graphene [6] using a carbon-oxygen van der Waals potential [11], or properties of biomolecules [7] using the GROMOS53A6 force field [12] (a carbon-oxygen Lennard-Jones interaction) or a Lennard-Jones potential using varying parameters in [5].In the FFMD simulations carried out in this study, we used three different water-carbon interaction potentials, with two 2-body potentials, one adjusted on DFT which was used for the data reported in the main text [13] and one on [14], and a carbon wall potential adjusted on water adsorption from quantum Monte Carlo data [15].Figure S5 clearly shows that our simulations gives similar results as reported in literature, with all FFMD simulations yielding a similar decrease of ϵ ⊥ as a function of the water slab thickness, irrespective of the force field used.In particular, they all match for water under strong confinement, when the distance between the confining surface is less than 1 nm.This indicates that when the strength of the confinement becomes large, the exact details of the interaction have no impact.Indeed, under strong confinement, the orientation of the water molecule is constrained in the in-plane direction due to the small space available in the slit, thus the ability of the water dipole to reorient with the external electric field is strongly suppressed.For larger water slabs, the shift in the predicted ϵ ⊥ associated with different force fields remains substantially negligible.Calculated out-of-plane dielectric constant of water confined in graphene slits using different force-fields in literature, reported in the following references: Brandenburg et al. [16], Ruiz-Barragan et al. [13], Steele [14], Loche et al. [7], Zhang [5], Itoh & Sakuma [9], Varghese et al. [6].Symbols labeled with a star (*) correspond to simulations made in this study, while the rest of simulated data are reported here from the references.Grey symbols are the experimental data reported here from Ref. [1].
We also investigated the effect of the water density on the calculated dielectric constant of the water confined inside graphene slits, as the water density in the slits used in the experiment is not known.To this aim, we run FFMD simulations with different values of water density ranging from 0.95 to 1.4 g/cm 3 , corresponding to lateral pressures between -0.3 kbar and around 6 kbar inside the confined region.Results are plotted in Fig. S6, again against the experimental data, clearly showing that the change in water density range considered here has only a minor impact on the predicted values of ϵ ⊥ of confined water.The small variation in ϵ ⊥ observed here arises from the competition Out-of-plane dielectric constant experimental data 0.95g/cm 3  1.00g/cm 3  1.10g/cm 3  1.20g/cm 3  1.40g/cm 3 Figure S6: Impact of bulk water density on the dielectric constant of confined water.Calculated dielectric constant of water confined in graphene slits for different values of bulk water density.Grey symbols are the experimental data reported here from Ref. [1]. between two effects.On the one hand, increasing the water density increases the number of molecular dipoles per volume unit.If the response of individual dipole moments remains unchanged, the polarization response is then expected to increase.On the other, increasing the water density would increase the rigidity of the hydrogen bonding network, which in turn decreases the ability of the water dipoles to reorient with the external electric field.Our simulations show a small decrease in ϵ ⊥ with increasing the water density up to 1.1 g/cm 3 , suggesting the increased connectivity of the hydrogen bonding network is dominant at relatively low density values.For density larger than 1.1 g/cm 3 , the trend reverses, indicating a small increase in ϵ ⊥ , thus the increase of dipole density becomes the dominant effect.However, the variations obtained here clearly show that both effects are not significant, with a decrease in ϵ ⊥ of only 10% for a change in water density of 20%.

S4 Hydrogen bond orientation
For the analysis of the hydrogen bonding network inside the slit described in the main text, we used the hydrogen bond definition based on a distance criterion between the two oxygen atoms r OO and an angle criterion on the θ HO d Oa angle, where O d is the oxygen of the donor molecule and O a is the oxygen of the acceptor [17].The two molecules form a hydrogen bond if r OO ≤ 3.5 Å, and θ HO d Oa ≤ 30 • .The direction of the hydrogen bond is then defined as the direction of the − −− → O d O a vector.Figure S7 compares the H-bond distributions obtained from NNP simulations, reported in Fig. 3 in the main text, with the distributions obtained from our FFMD simulations in the first two interfacial water layers near the graphene interface.From the comparison between the two data set, it can be seen that the obtained orientation is similar, but the NNP simulations yield a larger proportion of out-of-plane hydrogen bonds and in a smaller angular window in the first interfacial layer (L1).This indicates that a more rigid hydrogen-bonding network is predicted by NNP simulations as opposed to FFMD simulations.
Table S3 compares the average number of hydrogen bonds for the interfacial layers of water confined by carbon and hBN from NNP simulations.As observed already in FFMD simulations [18], the average number of hydrogen bonds is a non-monotonic function of the slit size, reaching a maximum for the bilayer case.This suggests a high structural order of water molecules in the bilayer.We note that the average number of hydrogen bonds for the interfacial layer L1 is found to be lower by 10% compared to what was reported for FFMD simulations [18].We also observed a significant drop in the number of hydrogen bonds when going from bilayer to monolayer, even larger than what observed in FFMD simulations [18].
Table S4 compares the average dipole moment of the water molecules in the interfacial layers confined by carbon and hBN from NNP simulations.We note that these results point towards a correlation between the average number of hydrogen bonds per molecule, and thus the local structure of water and the average dipole moment as seen for bulk water [19], which was expected from weakening of covalent bonds upon hydrogen-bond formation.
Figure S1: Computation of the dipole moment.Correlation between the dipole moment of the simulation cell computed from the Berry phase, the Wannier centers and the charge density.

Figure S2 :
Figure S2: Calculation of water's dielectric constant profiles.Snapshots (left) and corresponding dielectric profiles (right) for the total system (a), the system without water (b) and the profile (c) obtained after subtracting the polarization response of the graphene layer (left) and hBN layer (right) to the total system.

Figure S3 :
FigureS3: Convergence test of water's dielectric properties in AIMD simulations.Convergence of the cell dipole moment calculated using the Berry phase for an electric displacement of -1.0 V.nm −1 (a) with the corresponding running average of the total dipole moment for the M [1.935 nm] system and (b) running value of the dielectric constant of the interfacial water layers L1 and L2 for different water slab thicknesses.

Figure S4 :
Figure S4: Impact of the dividing surface.(a) Examples of positions of the dividing surface (dashed vertical lines) at the carbon atom side and calculated charge density profile using DFT calculations (solid line).A snapshot of the system is provided as a guidance.(b) Calculated inverse dielectric constant of water for various dividing surfaces shown in (a) using DFT and FFMD calculations.

Figure S5 :
FigureS5: Impact of the force-field choice.Calculated out-of-plane dielectric constant of water confined in graphene slits using different force-fields in literature, reported in the following references: Brandenburg et al.[16], Ruiz-Barragan et al.[13], Steele[14], Loche et al.[7], Zhang[5], Itoh & Sakuma[9], Varghese et al.[6].Symbols labeled with a star (*) correspond to simulations made in this study, while the rest of simulated data are reported here from the references.Grey symbols are the experimental data reported here from Ref.[1].

Figure S7 :
FigureS7: Orientation of hydrogen bonds in the first two interfacial layers.Orientation distribution with respect to the normal direction of the hydrogen bonds for water (a) in the first and (b) the second interfacial layer (L1 and L2) near the graphene surface, normalized by the isotropic distribution, calculated using NNP (solid lines) and FFMD (dashed lines) simulations for various slab thicknesses.The direction corresponds to the donor oxygen-acceptor oxygen.

Table S3 :
Average number of water hydrogen bonds in the first two interfacial layers.Number of the hydrogen bonds of the interfacial water layers molecules, L1 and L2, reported in Fig.1bin the main text, for both graphene (first column) and hBN slits (second column) and various water slab thickness using AIMD simulations.

Table S4 :
Average molecular dipole moment of water in the first two interfacial layers.Calculated molecular dipole of the water molecules in the interfacial layers L1 and L2, reported in Fig.1bin the main text, for both graphene (first column) and hBN slits (second column) and various water slab thickness using AIMD simulations.