# Fileset [s41467-025-58608-6.pdf](https://mdr.nims.go.jp/filesets/a574057a-c009-4d02-b52e-a0cf3e23acd8/download) ## Creator Ziwei Wang, Anupam Bhattacharya, Mehmet Yagmurcukardes, Vasyl Kravets, Pablo Díaz-Núñez, Ciaran Mullan, Ivan Timokhin, [Takashi Taniguchi](https://orcid.org/0000-0002-1467-3105), [Kenji Watanabe](https://orcid.org/0000-0003-3701-8119), Alexander N. Grigorenko, Francois Peeters, Kostya S. Novoselov, Qian Yang, Artem Mishchenko ## Rights [Creative Commons BY Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0/) ## Other metadata [Quantifying hydrogen bonding using electrically tunable nanoconfined water](https://mdr.nims.go.jp/datasets/0f346be1-84ab-4ae0-ac16-c16890716673) ## Fulltext Quantifying hydrogen bonding using electrically tunable nanoconfined waterArticle https://doi.org/10.1038/s41467-025-58608-6Quantifying hydrogen bonding usingelectrically tunable nanoconfined waterZiwei Wang 1,2 , Anupam Bhattacharya1, Mehmet Yagmurcukardes3,Vasyl Kravets1, Pablo Díaz-Núñez 1,2, Ciaran Mullan1, Ivan Timokhin1,2,Takashi Taniguchi 4, Kenji Watanabe 4, Alexander N. Grigorenko 1,Francois Peeters 5, Kostya S. Novoselov 1,2,6, Qian Yang 1,2 &Artem Mishchenko 1,2Hydrogen bonding plays a crucial role in biology and technology, yet itremains poorly understood and quantified despite its fundamental impor-tance. Traditional models, which describe hydrogen bonds as electrostaticinteractions between electropositive hydrogen and electronegative acceptors,fail to quantitatively capture bond strength, directionality, or cooperativity,and cannot predict the properties of complex hydrogen-bonded materials.Here, we introduce a concept of hydrogen bonds as elastic dipoles in anelectric field, which captures awide range of hydrogen bonding phenomena invarious water systems. Using gypsum, a hydrogen bond heterostructure withtwo-dimensional structural crystalline water, we calibrate the hydrogen bondstrength through an externally applied electric field. We show that ourapproach quantifies the strength of hydrogen bonds directly from spectro-scopic measurements and reproduces a wide range of key properties of con-fined water reported in the literature. Using only the stretching vibrationfrequency of confined water, we can predict hydrogen bond strength, localelectric field, O-H bond length, and dipole moment. Our work also introduceshydrogen bond heterostructures – a class of electrically and chemically tun-able materials that offer stronger, more directional bonding compared to vander Waals heterostructures, with potential applications in areas such as cata-lysis, separation, and energy storage.Hydrogen bonds (HBs) are unique intermolecular interactions that arestronger and more directional than weak and isotropic van der Waalsforces, yet weaker and less directional than covalent bonds1–4. HBs aretypically defined as attractive interactions between a hydrogen atom(H) and an acceptor group (A) in a system represented as D-H···A,where H is covalently bonded to an electronegative donor group (D),while H carries a partial positive charge δ+5,6. A key feature of HBs istheir directionality, often with a D-H···A angle close to 180°. In waterand similar systems, where distance between D and A is typicallygreater than 2.5 Å, the electrostatic attraction betweenH (δ+) andA (δ−)is believed to be the primary factor contributing to the strength ofHBs7–11.HBs could significantly alter the properties of residing system. Forinstance, the formation of HBs in water has a pronounced impact onReceived: 7 June 2024Accepted: 24 March 2025Check for updates1Department of Physics and Astronomy, University of Manchester, Manchester, UK. 2National Graphene Institute, University of Manchester, Manchester, UK.3Department of Photonics, Izmir Institute of Technology, Izmir, Turkey. 4National Institute for Materials Science, Tsukuba, Japan. 5Department Physics,University of Antwerp, Antwerpen, Belgium. 6Institute for Functional Intelligent Materials, National University of Singapore, Singapore, Singapore.e-mail: ziwei.wang@manchester.ac.uk; qian.yang@manchester.ac.uk; artem.mishchenko@manchester.ac.ukNature Communications | (2025) 16:3447 11234567890():,;1234567890():,;http://orcid.org/0000-0003-4062-7978http://orcid.org/0000-0003-4062-7978http://orcid.org/0000-0003-4062-7978http://orcid.org/0000-0003-4062-7978http://orcid.org/0000-0003-4062-7978http://orcid.org/0000-0003-1511-6252http://orcid.org/0000-0003-1511-6252http://orcid.org/0000-0003-1511-6252http://orcid.org/0000-0003-1511-6252http://orcid.org/0000-0003-1511-6252http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0002-1467-3105http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0003-3701-8119http://orcid.org/0000-0002-4109-2672http://orcid.org/0000-0002-4109-2672http://orcid.org/0000-0002-4109-2672http://orcid.org/0000-0002-4109-2672http://orcid.org/0000-0002-4109-2672http://orcid.org/0000-0003-3507-8951http://orcid.org/0000-0003-3507-8951http://orcid.org/0000-0003-3507-8951http://orcid.org/0000-0003-3507-8951http://orcid.org/0000-0003-3507-8951http://orcid.org/0000-0003-4972-5371http://orcid.org/0000-0003-4972-5371http://orcid.org/0000-0003-4972-5371http://orcid.org/0000-0003-4972-5371http://orcid.org/0000-0003-4972-5371http://orcid.org/0000-0002-6203-7867http://orcid.org/0000-0002-6203-7867http://orcid.org/0000-0002-6203-7867http://orcid.org/0000-0002-6203-7867http://orcid.org/0000-0002-6203-7867http://orcid.org/0000-0002-0427-5664http://orcid.org/0000-0002-0427-5664http://orcid.org/0000-0002-0427-5664http://orcid.org/0000-0002-0427-5664http://orcid.org/0000-0002-0427-5664http://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-58608-6&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-58608-6&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-58608-6&domain=pdfhttp://crossmark.crossref.org/dialog/?doi=10.1038/s41467-025-58608-6&domain=pdfmailto:ziwei.wang@manchester.ac.ukmailto:qian.yang@manchester.ac.ukmailto:artem.mishchenko@manchester.ac.ukwww.nature.com/naturecommunicationsthe physical properties of its condensed phases—ice. Often, hydrogen-bonded systems formextensive networks ofHBs that extend over longranges and in various directions, as seen in liquid water, ices, hydro-gels, proteins, and DNAmolecules. The complexity of these networks,along with the presence of interfaces or external electric fields inrealistic heterogeneous systems, makes it difficult to analyze themusing classic definition of HBs or to predict the properties of suchsystems12. To address these challenges, the electrostatic model of HBs,initially proposed decades ago that based on point chargeapproximation13, has been refined with various treatments, such aselectron distribution calculations14 and orbital hybridization theory15,leading to a more sophisticated description of HBs. Despite advance-ments in modeling and the introduction of spectroscopic methods,which provide direct information on molecular vibrations ofhydrogen-bonded species16, a practical and efficient experimentalapproach to quantitatively predict the properties of HBs is still lacking.Here we propose an approach that conceptualizes an HB as anelectric dipolemomentp of the donor-hydrogen (D-H) pair interactingwith the electric field EHB induced by the acceptor A (Fig. 1a). Thestrength of EHB at the dipole’s position depends on the electro-negativity of the acceptor A, and the distance between EHB and p.Considering that the typical distance between A and the D-H dipole inan HB is oftenmore than twice the length of the dipole, we assume EHBto be uniform at the dipole’s position. Also, the EHB is the effectiveelectric field from the acceptor, where the influence on it from p hasbeen embedded in ourmodel. The dipolemoment (p = qd) arises fromthe separation of the charge q along the D-H bond direction dD→H,while the field EHB generated by the acceptor nearby, A, provideselectric potential energy UHB = −p·EHB. This potential energy, whichdepends on themagnitude of p, the strength of EHB at dipole position,and cosine of the angle between them, represents the electrostaticinteraction energy of the hydrogen bonding system. It is maximizedwhen the electric field and the dipole are aligned. The field EHB stret-ches the D-H dipole, effectively lowering its spring constant from thatof anunbounded freeD-Hbond.When subject to external electric field(Eext), the projected component of Eext along the HB direction canfurther alter the spring constant and dipole moment of D-H dipole,thus changing the overall electrostatic interactions of the hydrogenbonding system. This offers amoreuniversal andprecisely quantifiabledipole-in-E-field alternative to the abstract chemical bond definition ofHBs. Ourmodel reflects both the dominant electrostatic nature of HBs(electric potential energy UHB as a quantitative measure of HBstrength) and their directionality (the angle between p and EHB,Fig. 1b). Moreover, by treating individual HBs as electric fields exertingon separate dipoles, our method allows the overall properties of asystem with complex HB networks to be predicted based on thesuperposition of dipole moments and electric fields. In this work, weintroduce the dipole-in-E-field approximation for HBs and use� ���d- d+ d-���� �������������� ���� ����������� ���� ������ � ��������������������������� Fig. 1 | Hydrogen bond model and gypsum system. a Hydrogen bond models.The upper diagram represents conventional model. Red and green atoms labeledwith D and A are negatively charged (δ−) donor and acceptor, respectively.Hydrogen (H) atom carries partial positive charge (δ+). The lower diagram repre-sents our dipole-in-E-field approach, where the HB is conceptualized as an elec-trostatic interaction of a dipole oscillator, with spring constant k and dipolemoment p (blue arrow), in the local field EHB (yellow arrow) of a negatively chargedacceptor. b Schematics of the HB in (a) when subjected to an external electric field(Eext) along the out-of-plane direction. The orientation of the dipole is described byangles ϑ and ϕ, which is the angle of D-H bond with respect to in-plane crystallineaxis and out-of-plane confining axis. c Unit cell of gypsum (in the box) and theextended confinedwater (oxygen in red, hydrogen in white) plane in the same two-dimensional molecular sheet. The crystal structure is simulated by our densityfunctional theory study. d Zoomed-in sketch of (c) showing two types of O-Hdipoles with different H atoms, intralayer O-HA and interlayer O-HB. e Ramanspectrumofwater stretchingmodes in an exfoliated thick gypsumcrystal, assignedto the stretching of O-HA and O-HB dipoles. Schematics of stretching vibrations areshownnext to the respective peaks as insets. The top inset is an opticalmicrographof exfoliated thick gypsum flake on a CaF2 substrate. The polarization of laser isdescribed by ϕE, the angle between laser polarization (green arrow) and gypsumcrystalline c-axis (dashed line). Scale bar 30 μm. The spectrum was taken atϕE = 90°. f Angle-resolved Raman intensity of O-HA (blue squares) and O-HBstretching (green circles) stretching modes. Solid lines the sinusoidal fittings usedto obtainϕ for O-HA andO-HB (ϕ defined in (b)) of water molecules in gypsum. Theprojection of a watermolecule on the gypsum a-c plane is shown in the center as aninset. Source data are provided as a Source Data file.Article https://doi.org/10.1038/s41467-025-58608-6Nature Communications | (2025) 16:3447 2www.nature.com/naturecommunicationsnaturally occurring hydrogen-bonded heterostructure gypsum toquantitatively describe HBs in nanoconfined water.ResultsDipole-in-E-field model of hydrogen bondingConventionally, the strength of HBs is characterized as a shift in thestretching vibration frequency of the D-H covalent bond: the strongerthe HB, the more red-shifted the frequency (ωD-H). Here we considerthe HB as D-H dipole in the presence of a local electric field, with thered shift in frequency depending on the local field produced by theacceptor (EHB). In a harmonic first-order approximation, the totalpotential energy (Uosc) and electric potential energy (UHB) of D-Hdipole oscillator in E-field are given by:Uosc =12kð xj j � d�� ��Þ2 +UHB ð1ÞUHB = � p � E ð2ÞHere, k is the force constant of the D-H bond, d and x are vectorspointing from the donor (at rest) to the equilibrium and instantaneouspositions of hydrogen, respectively. The E field is the effective electricfield at the dipole position (see Methods, ‘Full form of the dipole-in-E-field approximation’), which changes the equilibrium bond length andbond alignment elastically. In the absence of rotation, the problem issimplified to one dimension, with the vectors d, x, and E reduced totheir colinear components along D-H bond (seeMethods, ‘Full formofthe dipole-in-E-field approximation’). Using Hooke’s law, d(E) = d(0) +F(E)/k(E), we can express both the electrostatic force F(E) and the forceconstant k(E) as derivatives of the potential energy around equilibriumas F Eð Þ= � ∂UHB∂d and k Eð Þ= � ∂2Uosc∂d2 , leading, in a linear approximation,to the E-field dependent equilibrium bond length d(E) and forceconstant k(E)17:d Eð Þ=d 0ð Þ+ Ek Eð Þ∂p∂dð3Þk Eð Þ= kð0Þ � ∂2p∂d2 E ð4ÞHere, d(0) and k(0) are the equilibrium bond length and the forceconstant of a free D-H bond in the absence of the HB, respectively. Inthe harmonic approximation, the force constant k(E) of D-H bond canbe determined experimentally (using vibrational spectroscopy) fromthe stretching vibration frequencyωD-H of the D-H bond, k Eð Þ=μω2D�H,where μ is the reducedmass of the D-H system. Equations 3 and 4 givethe first and second derivatives of the D-H dipole moments withrespect to distance, enabling us to calculate the dipole moment of thehydrogen-bonded D-H system from Taylor expansion, and to deriveparameters related to hydrogen-bonding confinement such as HBenergy and dielectric properties.Experimental realization of the model in gypsumTo obtain the coefficient ∂2p/∂d2 in the Eq. 4, we performed vibra-tional spectroscopy experiments onwater HBs in a solid-state systemwhere the rotational degree of freedom of water molecules is sup-pressed. We chose gypsum (CaSO4·2H2O) as our model system.Gypsum is one of the most abundant minerals on Earth and consistsof alternating water bilayers and ionic CaSO4 sheets18,19. In gypsum,water molecules are confined and orientated between CaSO4 wallsvia HBs, forming a two-dimensional (2D) water crystal that extendsalong the gypsum a–c plane (Fig. 1c). We selected gypsum for threemain reasons. First, the low dimensional hydrogen bonded water ingypsum offers a favorable platform to study 2D confined water – ahot topic by its nature, and simplifies relevant theoretical analysis.Second, the stretching vibrations of the two O-H dipoles in thegypsum water molecules are sufficiently decoupled. In gypsum,water molecules are arranged in four orientations with two distincttypes of HBs, resulting in two groups of O-H dipoles in differentcrystallographic environments19,20: O-HA, bonded by a strongerintralayer HB, and O-HB, with a weaker interlayer HB nearly parallel tothe b-axis (Fig. 1d). Unlike the liquid water where the O atoms inwater serve as both hydrogen donor and acceptor within small waterclusters21, the distance between O atom in one water molecule and Hatom in another exceeds 2.8 Å, much greater than the typical lengthof H…O in conventional HB, suggesting minimal interactionsbetween water molecules. Therefore in gypsum, the O atoms in onewater molecule do not serve as HB acceptor for another, indicatingthe absence of intermolecular HBs within the 2D water layer. Thewater molecules in gypsum are bonded only to the CaSO4 molecularwalls. This has been confirmed by our DFT simulation. Additionally,we neglected the interaction between Ca2+ ions and O atoms in waterdue to the considerable distance separating them (~2.3 Å). Figure 1eshows the Raman spectrum of the H2O stretching vibrations in a bulkgypsum crystal, with two peaks centered at 3405 cm−1 and 3490 cm−1attributed to the stretching vibrations of O-HA and O-HB bonds,respectively (see details in Methods, ‘Decoupling of the stretchingfrequencies of O-H bonds in crystalline water of gypsum’). The polarplot in Fig. 1f shows the angle-resolved Raman intensity of O-HA andO-HB modes with respect to the a-c plane of gypsum. The angles ofmaximal intensity are ϕE = 18° for O-HA and ϕE = 90° for O-HB, with adifference of 72°. In subsequent experiments, we used these polar-izations to obtain maximum O-H stretching signal.Third, gypsum can be easily exfoliated into few-layer sheets,which are stable in ambient conditions. This enables us to apply vander Waals technology22 to construct gypsum heterostructures, whichallows us to apply external electric field Eext while simultaneouslyprobing the vibrational response of confined water. To this end, wesandwiched gypsum crystals (50–150nm thickness) between hex-agonal boron nitride (hBN) and few-layer graphene (FLG) to formheterostructures (Fig. 2a, b), resembling a parallel plate capacitor. Thetop and bottom FLGs allow us to apply Eext as a perturbation to thelocal field EHB of the acceptor A (in the case of gypsum, A is the oxygenatom in sulfate group, see Fig. 1b).We thenmeasuredRamanspectraofthese twoO-Hstretching vibrationswith the Eextfield applied along theout-of-plane b-axis.Electrically tunable stretching vibration modesThe color maps in Fig. 2c, d show the evolution of the Raman spectrafor O-HA and O-HB stretching vibrations as a function of Eext at 80K,respectively. Both maps illustrate the impact of Eext on the O-Hstretching vibration frequency, with O-HA showing pronounced peaksplitting and O-HB showing peak shifting. These observed changes areattributed to alterations of the total electric field (Etot = EHB + Eext)acting on different types of O-H bonds in gypsum. The inset of Fig. 2edemonstrates how Etot differs for O-H dipoles that point towardsopposite directions. We observed similar behavior for the O-HAstretching in Fourier transform infrared (FTIR) spectroscopy mea-surements at low temperature (Supplementary Fig. 1). The observedeffects from Eext are reversible within our measurement range of±0.5 V nm−1, as the spectra revert to their original state when Eext wasremoved (Supplementary Fig. 2). This indicates that the field-inducedchanges in the spectra are elastic effects rather than permanentstructural changes to gypsum. We also ruled out the possibility ofspectrum evolution due to electrostatic pressure effects, as the pres-sure induced by Eext is estimated to be ~0.01GPa at Eext = 0.5 V nm−1.According to reported values23, such pressure would cause shifts ofonly 0.01 and 0.2 cm−1 for O-HA and O-HB, respectively, which arebelow our detection limit.Article https://doi.org/10.1038/s41467-025-58608-6Nature Communications | (2025) 16:3447 3www.nature.com/naturecommunicationsFigure 2e summarizes the dependence of the Raman peak posi-tions of the O-HA and O-HB modes on Eext. The procedures forextracting peak positions are detailed in Methods, ‘Data analysis ofRaman spectra.’ O-HA stretching mode peak position is temperatureinsensitive, since the Eext-induced splitting remained unchangedacross the measured temperatures (10, 80, and 300K). For O-HBstretching mode, however, a red-shift from 3497 to 3481 cm−1 wasobserved at zero Eextwhen temperature decreased from300K to 80K.The FTIR peak splitting of the O-HA stretching is summarized in Sup-plementary Fig. 1. The difference between Raman and IR peak posi-tions at zero Eext arises from the frequency mismatch of the Raman-active mode (Ag) and IR-active mode (Au), as confirmed by our DFTcalculations (see Methods, ‘Decoupling of the stretching frequenciesof O-H bonds in crystalline water of gypsum’ and SupplementaryFig. 3). We also investigated the H-O-H bending mode under Eext usingscattering-type scanning near-field optical microscopy, as the bendingmode provides direct access to the HB characteristics24. However, noEext-induced peak splitting or shifting was detected, consistent withour DFT results (see Methods, ‘Water bending mode in gypsum’ andSupplementary Fig. 4 for details).Fitting dipole-in-E-field model with experimental dataNow, we recall the model at the beginning of the paper to rationalizeour experimental results. First, we convert the Eext-dependentstretching mode peak positions (ωD-H) into force constants,k(E) =μω2D�H, to obtain k(E) as a function of Eext. The total electricfield acting on an O-H dipole (denoted as Etot), is the sum of theinternal E-field (EHB) and the external electric field projected on theO-H dipole (Eext cos ϑ):Etot = EHB + Eext cosθ ð5ÞHere, θ is the angle between and the gypsum b-axis (Fig. 1b). For thedipoles O-HA and O-HB in gypsum, θ is 69.60° and 9.37°, respectively,as obtained from neutron diffraction data20. For a free watermolecule,k(0) = 762.0 ±0.3 Nm−1 25. Fromhere, we construct a plot of k(E) vs Etot,Fig. 3a. Further details of data processing and error analysis are inMethods, ‘Fitting of k vs E relation in Fig. 3a and local fields at dipolesO-HA and O-HB.’ This analysis allows us to quantify the internal electricfield EHB forO-HA andO-HB as 5.33 ± 0.63 V nm−1 and 3.82 ± 0.45 V nm−1,respectively. From Eq. 4, we obtain the slope of k(E) vs Etot curve inFig. 3a, the second derivative of the dipole moment along the O-Hbond, ∂2p∂d2 = (2.2 ± 0.3) × 10−8Cm−1.Similarly, we further obtain the first derivative of the dipolemoment, ∂p∂d, using Eq. 3. Todo this,we rely onexisting literaturedataofseveral solid-state systems containing crystalline water, such asproton-ordered and proton-disordered ices, and solid hydrates. Forfree water molecule, we used the equilibrium O-H bond lengthobtained from fitting databases of experimental rotational energylevels26,27. For neutron diffraction data we applied thermal corrections,more details are in Methods, ‘Thermal corrections for neutron struc-tural data.’ From here, we establish the d(E) vs E/k(E) relation, Fig. 3b,Fig. 2 | Electric field effects on O-H dipole stretching vibration. a Schematic ofRaman measurements with external electric field, Eext, applied to a tgr/hBN/gyp-sum/hBN/bgr heterostructure. b Optical micrograph of one of our gypsum het-erostructure devices. tgr and bgr refer to top and bottom FLG electrodes. Scale baris 20μm.Colormaps of Raman spectra forO-HA (c) andO-HB (d) at 80K under Eext.Color from dark blue to yellow indicates increasing Raman intensity. e Raman peakpositions of the two O-H stretching modes at different temperatures under Eext.Datasets with the same temperature but different colors were collected from dif-ferent samples. Six datasets were acquired from six separate samples. The insetshows the change in O-H dipoles and Etot under Eext, the effect is exaggerated forclarity. Dashed circlesmark positions of H atoms in the absence of Eext. Source dataare provided as a Source Data file.Article https://doi.org/10.1038/s41467-025-58608-6Nature Communications | (2025) 16:3447 4www.nature.com/naturecommunicationsusing systems for which literatures contain both neutron diffractiondata to obtain O-H bond length d, and Raman spectra to get kE and Efrom O-H stretching vibration frequency. Linear fitting of Fig. 3b givesthe first derivative of dipole moment, ∂p∂d = (7.7 ± 1.2) × 10−19C (Eq. 3),and the fitted O-H bond length of a free water molecule asd(0) = 0.9570Å (literature sources, details of data analysis, as well asfitting details are in Methods, ‘Data analysis for bond length vs E/krelation in Fig. 3b’). Given that both data from solid hydrates and dif-ferent phases of ices have been used in the fitting in Fig. 3b, ourmodelis designed to be applicable to both solid hydrate and ice systems.We then reconstruct a universal expression for O-H dipolemoment p(d) using the obtained derivatives in a Taylor expansion:p dð Þ=X1n =01n!∂np∂dn d � d 0ð Þð Þn ð6ÞBy using the first two derivatives (∂2p∂d2 and ∂p∂d) obtained fromFig. 3a, b, and the knownO-H dipolemoment for a freewatermoleculep(d(0)) = 1.5140±0.0004D26–31, weplot theO-HdipolemomentpO-H asa function ofO-Hbond length in Fig. 3c. The calculateddipolemomentpH2O = 1.88 ±0.04D for free water molecule agrees well with thereported value pH2O= 1.855 ± 0.001 D29. For liquid water, our modelpredicts pH2O = 2.66 ±0.13 D, aligning well with the reported values of2.9 ± 0.6D32withinuncertainty range. For ice Ih, ourmodel estimated adipole moment of 3.39 ± 0.24D for a proton-order ice Ih structure,which is higher than the dipole moment reported for proton-disordered ice Ih, 3.09 ±0.03D, derived from the multipolemoments expansion of water molecule33. The overestimated dipolemoment may result from the absence of defects that are typicallypresent in ice Ih. Thesedefects facilitate the formationofH3O+ andOH−pairs and enable proton hopping in ice Ih, which canweaken theHBs34.To exclude the effect from defects, we compared our results withthose from the proton-ordered counterpart of ice Ih, commonlyknown as ice XI. The dipole moment of 3.39 ±0.24D estimated by ourmodel aligns closely with the 3.3 D reported for water molecules in iceXI, given by DFT simulations35. All estimated dipole moments togetherFig. 3 | Hydrogen bonding and dielectric properties of water. a Force constant kvs Etot along O-H bond. Symbols represent experimental results from this workusing gypsumatdifferent temperatures. The linearfit using Eq. 4 (red line) gives thesecond derivative of the dipole moment along the O-H bond, ∂2p/∂d2. b O-H bondlength d versus Etot/k for a range of solid-state water systems from literature(references are inMethods, ‘Data analysis for bond length vs E/k relation in Fig. 3b’).Error bars for dwere provided in the references; error bars for Etot/kwere obtainedfrom the calculations using our fitted ∂2p/∂d2 from (a). Linear fit using Eq. 3 (purpleline) gives the first derivative of the dipole moment, ∂p/∂d. c Estimated dipolemoment along O-H bond, pOH(d), constructed using Eq. 6 (gray line), plottedalongside with different water systems from the literature (for details and refer-ences see Methods, ‘Dipole moments vs O-H bond lengths in Fig. 3c’). d Effectivepolarizability along the confining direction (αc) of a single confined O-H dipoleagainst theO-Hbond length d and the angle ϑbetween the dipole and the confiningdirection. α is the polarizability of O-H along dipole direction (ϑ =0°) given by Eq. 7.Red, blue, green, and gray dots represent the O-H dipoles in proton-ordered ice XI,II, IX, and gypsum. The arrows point to the structures of ice XI, II, and IX, with thegeometrically non-degenerate water molecules highlighted. Source data are pro-vided as a Source Data file.Article https://doi.org/10.1038/s41467-025-58608-6Nature Communications | (2025) 16:3447 5www.nature.com/naturecommunicationswith reference values are summarized in Supplementary Table 2. Fordetails about dipole moment calculations, see Methods, ‘Dipolemoments vs O-H bond lengths in Fig. 3c.’Dipole-in-E-field model for hydrogen bond energyFurther,we estimate theHBenergyUHB (Eq. 2) using thederiveddipolemoment and electric field values, given by Eqs. 6 and 4, respectively.For example, Ice Ih has a Raman shift of 3279 cm−1 and a bond length of1.014 Å (see Methods, ‘Data analysis for bond length vs E/k relation inFig. 3b’). Liquid water has a bond length of 0.984Å (see Methods,‘Dipole moments vs O-H bond lengths in Fig. 3c’) and a broad O-Hstretching band due to the extensive HB network, ranging from3200 cm−1 to 3650cm−1 36,37. Using our approach, we obtainUHB per HBas 10.5 ± 1.2 kcalmol−1 (455 ± 52meV) for ice Ih and 5.1 ± 0.5 kcalmol−1(221 ± 22meV) for liquid water. Within the uncertainty range, our UHBfor liquidwater is consistentwith 5.26 kcalmol−1, the reference value ofthe liquid water HB strength estimated from the enthalpy of vapor-ization of water at 25 °C38. However, our estimation for ice Ih is higherthan the reported lattice energy, 7.04 kcalmol−1 (305meV)39,40. Thisdifference could be due to an overestimated dipole moment calcu-lated on a proton-ordered ice Ih structure where H atoms have fixedpositions, as discussed above. In the real ice Ih, protons are disorderedand can hop between different crystalline sites, which tends to weakentheHB strength. Suchweakening effect is reflected in ourmodel by theoverestimated UHB for ice Ih. All estimated UHB and the comparisonwith reported values are summarized in Supplementary Table 2.Predicting the dielectric properties of HB systemsBeyond HB characterization such as dipole moment and HB energy,our approach could also predict other properties of the HB system,such as polarizability and dielectric constant. Below we show brieflythe prediction power of our approach. For example, polarizabilityalong the O-H dipole, α(d) = ∂p/∂E, can be evaluated using the rela-tionship d(E) from Eq. 3, as,α dð Þ= ∂p∂d+∂2p∂d2 d � d 0ð Þð Þ" #3=∂p∂dkð0Þ ð7ÞHere, ∂p∂d and∂2p∂d2 derived from our model, and d(0) and k(0) taken fromthe free water molecule (see Methods, ‘Fitting of k vs E relation inFig. 3a and local fields at dipole O-HA and O-HB’). In nanoconfined andcrystalline water systems, the rotational degrees of freedom of thedipole are restricted, meaning that no dipole reorientation under Eextwill contribute to the dielectric response. Therefore, Eq. 7 fullydescribes the polarizability of such systems. The effective polariz-ability along b-axis is given by αc = α cos2θ and is determined by both dand ϑ, as summarized in Fig. 3d, along with exemplar nanoconfinedwater systems. For system with known bonding geometry and bondlength d, the polarizability of each O-H dipole along any confining axiscan be evaluated via Fig. 3d. And the polarizability along the confiningdirection can be obtained by summing up αc for all dipoles in thesystem.From here, we further extract the static dielectric constant εr of asystem consisting of electrostatically confined water molecules (seealso Methods, ‘Polarizability and dielectric behavior of confinedwater’). In such a system, the polarizabilities of individual O-H dipolesadd up to give εr:εr = 1 +1ε0XiαciNi ð8Þwhere αci is the effective polarizability along confining direction andNiis the number density of the ith type of O-H dipoles. Using this relation,we obtain εr = 6.6 ± 0.7 for the confined 2D water layer in gypsum. Thelow εr of the 2D water sheets is consistent with the decreasing εr ofconfined liquid water along perpendicular direction, due to restrictedrotational freedom41,42. Our capacitance measurements of gypsumcrystal give εr = 3.68 (see Methods, ‘Polarizability and dielectricbehavior of confined water’ and Supplementary Fig. 5). Since gypsumis an hydrogen bond heterostructure (HBH) composed of alternatingCaSO4 and water dielectrics, we can evaluate εr = 2.8 ± 0.5 for the bareCaSO4 molecular sheets, whichmatches well with the reported εr of β-anhydrite, 2.7 ± 0.0243. The estimated εr for gypsum system issummarized in Supplementary Table 2.Our model also accurately replicates the dielectric constants forthree proton-ordered ices: ice XI (ordered ice Ih), ice II, and ice IX(ordered ice III). The polarizability ofO-Hbonds is given by αc = α cos2θandEq. 7, as shown in Fig. 3d (see alsoMethods, ‘Dielectric constants ofwater ice systems’). For bulk liquid water, using the dipole momentfrom Eq. 6 and considering rotational averaging, we obtain εr = 27 ± 2,which is much lower than experimental value of 79 (see Methods,‘Dielectric constant of liquid water’). This discrepancy arises becausethe rotations of watermolecules are not independent of each other, aneffect usually accounted for by introducing the dipole correlationfactor Gk44. To reproduce the experimentally obtained εr, Gk ≈ 3.15 isneeded in our case, falling in the range of calculated value of 2.72−3.70for bulk liquid water45. Ourmodel also cannot predict the anomalouslyhigh dielectric constant in proton-disordered ices, such as εr on theorder of 100 for ice Ih, themost common formof ice in nature46,47. Thisbehavior is usually attributed to the presence of defects in the form ofpartial rotations and proton jumps, which form polarized domains atlower temperatures where ice sometimes becomes ferroelectric48.Our approach also reproduces the anomalously low εr,⟂ of con-fined water, which arises from the restricted rotational degrees offreedom of liquid water under nanoconfinement49. Here, we considerthe dielectric response of confined water resulting only from electro-nic and atomic polarization, as described by Eq. 8. The dependence ofεr along the confining direction on the orientation ofwatermolecule isplotted in Supplementary Fig. 6, calculated using the density of con-fined water 3.39 g cm−3 from recent simulations50. We found that whenthe water molecular plane is nearly parallel to the confining plane(confining angle θ < 23°), we can reproduce the experimentallyobtained εr,⟂ = 2.149.Beyond conventional HB systemsOur dipole-in-E-field approximation can be generalized to other typesof electrostatic interactions beyond conventional HBs. Instead ofconventional chemical bonds, our model conceptualizes HBs as uni-versal electrostatic interaction between anO-Hdipole andelectricfieldfrom acceptor. Consequently, although the experimental data wascollected on solid hydrate like gypsum, our approach can be gen-eralized to any systemsbased on electrostatics interactions containingO-H dipole, including conventional hydrogen bonding systems likeliquid water and ice, and π-hydrogen bonding systems such as water-carbon interface. Once the vibrational frequency of water under con-finement is known, the strength of the electrostatic interactionbetween water and its surroundings can be evaluated. For example,interfaces between water and carbon materials often exhibit unex-pected phenomena such as low friction or long slip length, which arepromising for water-energy applications51. Such electrostatic interac-tions between water O-H bonds and the π orbitals of sp2-hybridizedcarbon canbe viewed asweakHBs, evidencedby red-shifts in thewaterstretching vibration frequency. Using our approach, we can evaluatethe interaction energy between water O-H dipoles and π orbitals atgraphitic surfaces.For instance, in the caseof singlewatermolecules trapped inside aC60 fullerene52, from the measured symmetric v1 and antisymmetric v3stretching modes at 3573 cm−1 and 3659.6 cm−1, we can calculateinteraction strength UHB = 1.5 ± 0.2 kcalmol−1 per bond (66 ± 9meV),Article https://doi.org/10.1038/s41467-025-58608-6Nature Communications | (2025) 16:3447 6www.nature.com/naturecommunicationscomparable to typical weak HBs. Similar values are obtained for waterconfined in single-walled carbon nanotubes53. Further evidence sup-porting our approach comes from recent sum frequency generation(SFG) measurements at the water-graphene interface54, where theauthors reported vibration frequency of νOH = 3616 cm−1 for the dan-glingO-Hdipole. Our approachpredicts theHBenergy in this case veryclose to that for water confined in C60. More conclusively, thestretching frequency of this dangling O-H bond shows an electric fieldtuning effect that matches perfectly with our model (see Methods,‘Water-carbon interface.’ and Supplementary Fig. 7). More details onevaluating water-carbon interactions see Methods, ‘Water-carboninterface.’DiscussionOur dipole-in-E-field approach reliably quantifies the characteristicsof confined water molecules in various complex environments,advancing our understanding of HB in biologically and technologi-cally relevant systems. Using coefficients parametrized from thiswork, we can quantitatively evaluate properties such as O-H bondlength, local electric field, dipole moment, HB strength, polariz-ability, and dielectric constant directly from spectroscopic data(Raman, FTIR, or SFG) for a wide range of water systems, includingcrystallohydrates, ices, and nanoconfined water. This approach canbe extended beyond confined water to systems such as organicmolecules, (bio)polymers, liquid phases, and electrochemical inter-faces, which are crucial for bio-medicine and technology. Testing thelimits of our approach in scenarios where quantum and many-bodyeffects dominate, such as blue-shifting strong HBs7 or cooperativeHB systems55, could provide new insights into the fundamental nat-ure of hydrogen bonding.Notably, our work represents the first study of hydrogen bondheterostructures (HBHs). While previous research has explored het-erostructures composed of hydrogen-bonded organic frameworks(HOFs) and graphene, highlighting how intralayer HBs within theHOF enhance interactions with graphene56, our study primarilyfocuses on the interlayer HB that binds layers along out-of-planedirection, which could emerge as a class of intelligent materialsallowing for facile assembly and fabrication, similar to van der Waalsmaterials. The highly anisotropic HBHs—which potentially can beprobed using dichroic IR measurement57—offer unique properties,including electrical and chemical tunability, stronger and moredirectional bonding compared to van der Waals forces, and thepotential for self-assembly and emergent phenomena. By com-plementing hydrogen-bonded organic frameworks, HBHs and HBtechnology could enhance applications in gas separation, catalysis,and (opto)electronic devices58–63. The combination of our dipole-in-E-field approach and the development of HBHs opens new avenues forthe rational design and precise control of the properties ofhydrogen-bonded systems.MethodsFull form of the dipole-in-E-field approximationThe 3D form of Hooke’s law can be written as F= � κ � X. Here X is thedisplacement vector in 3D space, and κ is a 3 × 3 strain stiffness tensor,where the diagonal and off-diagonal indices are the stiffness in direc-tions parallel and normal to each force components. The potentialenergy of a harmonic oscillator in 3D case is:Uosc =12XT � κ � X +UE , ð9ÞUHB = � p � E+ . . . ð10ÞSet x-axis along D-H bond, then we have electrostatic forceand sum potential stiffness: F(E) = −∇UHB and k Eð Þ=∇2Uosc. Thefield-dependent bond length and field-dependent spring constant arethen:d Eð Þ=d 0ð Þ+∇p � E 1k Eð Þ ð11Þk Eð Þ= kð0Þ � ∇2p � E ð12ÞIn this work, we consider only the parallel response along D-Hbond. Thus, among the nine indices of κ, only one diagonal term isnon-zero. The system then simplifies to the 1D problem described inthe main text.The electric field in our model is defined as the effective electricfield at the dipole position, generated by the hydrogen acceptor, A,with the modification from dipole-field interactions inherentlyaccounted for. As the dipole moves within acceptor’s field, the inter-action induces a charge redistribution in the acceptor, dynamicallyevolving the field and continuously affecting the HB strength. How-ever, these effects do not require separate evaluation. Using thedipole-in-field model, the HB energy is directly determined as theproduct of the effective dipolemoment and the effective electric field.Consequently, the modification of E due to dipole-field interactions isincorporated into our model by treating E in Eq. 2 as the effectiveelectric field at the dipole position, where the contribution of p hasbeen embedded. This effective E also corresponds to the Raman peakshift measured experimentally.Device fabricationGypsum films were isolated from bulk natural crystals onto SiO2(290 nm)/Si and CaF2 substrates (for spectroscopic measurements)and polydimethylsiloxane (PDMS) stamps (for heterostructure fabri-cation) using micromechanical exfoliation in cleanroom. Graphite/gypsum/graphite heterostructure is fabricated using all-dry viscoe-lastic stamping technique64: first, graphite flake (bottom electrode)was isolated from natural graphite crystal onto SiO2/Si substrate;gypsum flake and top graphite electrode were exfoliated onto PDMSstamp. The gypsum and top graphite flakeswere then transferred ontobottom graphite electrode successively, aligning and stacking targetflakes under the optical microscope using transfer rig.Gypsum is prone to dehydration in vacuum or when exposed toelevated temperature, therefore, conventional lithography and metaldeposition process cannot be used to fabricate electrical contacts togypsum. Instead, we used two thick graphite flakes with over 100μmlength as extended electrical contacts to the top and bottom FLGs ofthe heterostructure, to which gold wires were then attached usingconductive silver paint.For devices with hBN separation layers, the hBN flakes wereexfoliated and transferred using polypropylene carbonate (PPC) drytransfer technique65: the hBN flakes exfoliated on SiO2/Si wafer werepicked up using PPC coated PDMS stamp at 45 °C, and then alignedand transferred onto desired stacks at 60 °C under the optical micro-scope using transfer rig. All fabrication processes were carried out in acleanroom.Vibrational spectroscopy measurementsAngular-resolved polarized Raman spectroscopy was performed usingOxfordWITec alpha300R confocal Raman system,with 600 linesmm−1grating and a 532 nm wavelength laser with linear polarization. Laserpower was kept below 8mW. During measurements, the sample wasfixed on the sample stage, with incident laser polarization rotated by apolarizer at a step of 5°. The spectra were taken with parallel analyzerconfiguration at each incident orientation, so the collected light hadthe same polarization direction as incident. For each spectrum, theacquisition time was 10 s with 6 accumulations.Article https://doi.org/10.1038/s41467-025-58608-6Nature Communications | (2025) 16:3447 7www.nature.com/naturecommunicationsRaman measurements under electric field in ambient conditionswere carried out using Renishaw inVia confocal Raman system equip-ped with 2400 lines mm−1 grating and 514 nm linearly polarized laserwith power restricted below 8mW, focused through a × 100/0.90 NAobjective. The voltage was applied using a Keithley 2614Bsource meter.Peak splitting is not very pronounced even at 0.3 V nm−1 (e.g.,Supplementary Fig. 8 for O-HA stretching mode), but higher field isunreachable at ambient conditions due to field-induced degradationof FLG electrodes, even with hBN spacers. So, we performed Ramanmeasurement at T < 80K, where Eext of up to 0.5 V nm−1 can be appliedand distinct peak splitting is well resolved.Ramanmeasurements at cryogenic temperatureswere performedon WITec alpha300 R Raman system with 600 lines mm−1 grating and×100/0.85 NA objective using 532 nm linearly polarized laser withpower lower than 8mW.The samplewasplaced in anOxfordmicrostatsystem in vacuum with temperature controlled by a liquid N2 circula-tion loop formeasurements at 80K and liquidHe formeasurements at10K. The fluctuation of temperature was ±5 K. Due to the lightabsorption by the optical window on Oxford microstat system, theacquisition time was extended to between 30 and 60min, dependingon the signal to noise ratio on different devices, to increase the Ramanintensity. This significantly increased themeasurement duration, so nospatial mapping data was collected while the external electric fieldswere applied. To protect the samples from potential laser-induceddamage, the laser power for all Raman measurements was maintainedbelow 8mW.Fourier transform infrared (FTIR) spectra were recorded with aBruker VERTEX 80 spectrometer equipped with a Bruker HYPERION3000 FTIR microscope. For low temperature FTIR measurements at83 K, the sample was placed in a Linkam stage with temperature con-trolled by a liquid N2 circulation loop. All spectra were collected with4 cm−1 resolution and averaged over 512 scans.Nano-FTIR measurements were performed on an AFM basedscattering-type scanning near-field optical microscopy (s-SNOM) fromneaspec GmbH. The conventional Pt/Ir probe (ARROW-NCPt-50,Nanoworld) with a tapping frequency of 270 kHz was illuminated withap-polarizedmid-infraredbroadbanddifference frequencygeneration(DFG) laser (frequency range 1000–1900 cm–1, average power <1mW).The sample for s-SNOMmeasurements was FLG/gypsum/FLG stack onSiO2/Si substrate. During measurements, a dc voltage was applied tobottom graphite electrode, while the top graphite and AFM probewere electrically grounded, creating an electricfield alongout-of-planedirection. Obtained nano-FTIR spectra were collected at differentvoltage levels by averaging 30 interferograms (spectral resolution≈7 cm−1) with 2048 pixels and an integration time of 10ms/pixel. Allspectra are normalized to that of a Si substrate.Data analysis of Raman spectraFor Fig. 1f in themain text, we first performed background subtractionon raw spectra using 2nd polynomial with 256 sample points, followedby peak fitting using Voigt function. The peak position and peak widthof water stretching modes were constrained within 3400–3412 cm−1and 20–50 cm−1 for O-HA, and 3490– 3500 cm−1 and 20–70 cm−1 forO-HB. Raman spectra in Fig. 2 (main text) were obtained first throughbaseline fitting and subtraction using 4th polynomial with 256 samplepoints, followed by peak fitting using Voigt function. Initial peakpositions were set at the position with maximum intensity. The peakwidths were limited within 100 cm−1 and 90 cm−1 respectively.Decoupling of the stretching frequencies of O-H bonds in crys-talline water of gypsumIn gypsum, water O-H covalent bonds are stretched due to HB. Dif-ferent crystalline environment of the two O-H bonds results in a dis-torted asymmetric water molecule. Such asymmetry decouples andlocalizes the normal symmetric and antisymmetric stretching ontoindividual O-H bonds. Consequently, one O-H bond will have muchlarger vibration amplitude than the other in each vibration mode.We performed density functional theory (DFT) calculations usingQuantum ESPRESSO package, to elaborate on the vibrational responseof gypsum. Dispersion corrections that account for the van der Waalsforces were incorporated following the empirical D3 formulation byGrimme, yield 0.06% difference with experimental lattice parameters.Calculations of phonon frequencies were carried out using the finitedifference methodology implemented in Phonopy by displacing eachatom 0.001 Å along three principal axes. Generalized GradientApproximation (GGA) formulation by Perdew, Burke and Ernzerhofwas used tomodel exchange-correlation interactions among electronsin all the calculations. An energy cutoff of 65 Ry with a uniformMonkhorst-pack k-mesh of 5 × 5 × 5 is used for all the calculations onthe primitive unit cell. Convergence criteria used for all self-consistentfield calculations are 10−10Ry per formula unit.In gypsum crystal, water has 8 stretching modes in total, Sup-plementary Fig. 3. According to ourDFT calculations, they are groupedinto two branches, high- and low-frequency (near 3500 cm−1 and3400 cm−1), with vibration amplitudes mainly localized on O-HB andO-HAbonds, respectively. In eachof these twobranches, the four levelsarise due to the symmetry of gypsum crystal. In particular, four watermolecules in a unit cell are grouped into two inversion symmetricpairs. Each of the two water molecules in the pair vibrate collectivelyeither in-phase or out-of-phase, which decouples the modes. In turn,each in-phase and out-of-phase modes has symmetric and anti-symmetric intramolecular solutions, which leads to further modedecoupling.Infrared spectroscopy of water O-HA stretching modesin gypsumSimilar toRamandata presented in Fig. 2 inmain text, peak broadeningand splitting were also observed in low temperature FTIR measure-ments (Supplementary Fig. 1). The IR peak of O-HA mode blue-shiftedfrom 3401 cm−1 to 3412 cm−1 when temperature reduced from 300 to80K, possibly due to a volumetric contraction effect which shortensthe effective O-H bond length66. Similarly to Raman, the O-HAstretching splits into two peaks, corresponds to the twoO-HA Amodesin Supplementary Fig. 1b (AgA and AuA at zero field, SupplementaryFig. 3). The obtained electric field dependence of stretching peakposition in FTIR spectra is similar to that of our Raman results.Due to the large overlap between the two Eext inducedO-HApeaks,we instead chose positions of maximal peak intensity for spectra fit-ting. At high E-fields, both O-HA peaks show similar splitting withelectric field. However, at low E fields, the low-frequency peak is nearlyE-field independent. According to our DFT calculations, the two O-HApeaks correspond to two A modes of O-HA at 3406 and 3393 cm−1,which are only Raman-active (AgA mode) and IR-active (AuA mode) atzero field, respectively. They are not the samemode andhave differentfrequency. This accounts for the absence of AgA mode in IR and AuAmode in Raman at zero field, and the deviation of IR data from linearsplitting at low Eext-field (<0.065 V nm−1) in Supplementary Fig. 1c.Water bending mode in gypsumDue to extremely low Raman intensity, water bending mode and itsfield-dependence were investigated via nano-FTIR measurementsusing scattering-type scanning near-field opticalmicroscopy (sSNOM).Water in gypsum exhibits two bending modes located around 1620and 1680 cm−1, Supplementary Fig. 4. This agrees well with our DFTcalculations. The water bending mode also decouples into four sub-modes, because of four water molecules in a gypsum unit cell. Similarto the stretching vibration mode decoupling, inversion symmetrydivides these fourmodes into two pairs, with frequencies around 1620and 1680cm−1.Wedid not observe any E-field response for the bendingArticle https://doi.org/10.1038/s41467-025-58608-6Nature Communications | (2025) 16:3447 8www.nature.com/naturecommunicationsmodes, Supplementary Fig. 4. This is consistent with our DFT calcu-lations, which show that bending mode frequencies do not changeeven at a high electric field up to 1.8 V nm−1, because phonon disper-sion around Γ point was unaffected by external electric field.Fitting of k vs E relation in Fig. 3a and local fields at dipole O-HAand O-HBThe k vs E relation in Fig. 3a in the main text was obtained byaccounting for the applied Eext component with respect to O-Hdipoles and accommodating the crystal structure symmetry of gyp-sum, using data from Fig. 2e, main text. To this end, the Eext for eachdataset was scaled by cos ϑ to obtain the electric field componentalong each O-H dipole. Here ϑ (as defined in Fig. 1b, main text) is69.60° and 9.37° for O-HA andO-HB dipoles, respectively20. Then, halfof each dataset (those with positive slopes) were flipped horizontally,so that all the datasets having negative slope (red-shift). According toEq. 4 in main text, this slope corresponds to� ∂2p∂d2, where E is the totalelectric field at the dipole, equal to EHB + Eext. All datasets trans-formed by the above procedure were then fitted globally using Eq. 13with three fitting parameters—one is the shared slope ∂2p∂d2 for bothO-HA and O-HB data because ∂2p∂d2 is expected to be the same for bothdipoles, and two different EHB for O-HA and O-HB due to the differentlocal HB field:k Eð Þ= kð0Þ � ∂2p∂d2 ðEAðBÞHB + Eext cosθÞ ð13ÞTo obtain the intercept k(0) in Eq. 4 in main text (and Eq. 13), weresorted to the known data on water in gas phase from literature.Vibrational spectroscopy of freewatermolecules is complicated: threevibrational modes—symmetric v1 and antisymmetric v3 stretching andbending mode v2 together with numerous rotational modes result inhundreds of thousands of peaks of various intensities andlinewidths67–70. The most recent values of free H216O molecule’s vibra-tional transition frequencies from the W2020 database are 3657.053,3755.929, and 1594.746 cm−1, for v1, v3, and v2, respectively25. Thevibration energy of uncoupled O −H bond vO−H, estimated as theaverage of v1 and v3, is 3706.491 cm−1, in good agreement withvO−H = 3707.467 cm−1 for HD16O molecule, where decoupling is pro-videdby a largemismatch ineffectivemasses ofHandD71. This gives usthe harmonic spring constant k(0) ≈ 762.0 ±0.3 Nm−1 for uncoupledO −H bond in HB-free water. Linear fitting using Eq. 13 gives slope� ∂2p∂d2 = (2.23 ± 0.26) × 10−8Cm−1. The local EHB-fields at O-HA and O-HBdipoles are (5.33 ± 0.63) and (3.82 ± 0.45) V nm−1, respectively. Theuncertainty is given by the 2σ value of the linear fit (reduced χ2 = 2.3).Thermal corrections for neutron structural dataBond lengths from neutron diffraction are based on time-averagedneutron scattering intensity, giving an averaged interatomic separa-tion. The actual bond length is slightly longer, after applying thermalcorrections. The thermally correctedO-H bond lengths in Fig. 3b in themain text were obtained in different ways. For ice Ih, ice VIII, ice IX,Li2SO4·H2O, Ba(ClO3)2·H2O, andLiClO4·3H2O, their thermally correctedO-Hbond lengthswere provided in the respective references. For ice II,ice VI, ice VIII and ice XVII, the thermal corrections were evaluatedusing isotropic temperature factor, B. For CuSO4·5H2O andBaCl2·2H2O, thermal corrections were carried out using anisotropictemperature factors, βij. The temperature parameters of the systemsare given in corresponding references respectively. The approachesused to evaluated the thermal corrections were reported in ref. 72.Data analysis for bond length vs E/k relation in Fig. 3bThe various systems in Fig. 3b, main text contain ices and solidhydrates. The distance between H and metal ions (M) in most solidhydrates ranges from 1.5 to 3 Å. Here, we focused on pure hydrogenbonding interactions and analyzed the solid hydrates systems withH…M distance of 2–3 Å, ignoring materials such as BeSO4·4H2O,where the significantly shorter H…Be distance of 1.6 Å may affect thecharge distribution onH via electrostatic interactions fromBe cation.All O-H bond lengths are from thermally corrected neutron diffrac-tion data. Where the O-H length was measured on deuterated waterD2O, the thermal-corrected bond length is further upshifted by 3% toaccount for the bond difference between H2O and D2O73. For vibra-tional spectroscopy data on systems with multiple stretching peaks,we assignedpeaks toO-Hbondswith different bond lengths such thata shorter bond has a higher vibration frequency due to the steeperpotential profile. The assignment process for each confined-watersystem is discussed below and data is summarized in SupplementaryTable 1.For ice Ih, one Raman peak is observed at 3279 cm−1 for 3mol%HOD ice at 123 K, corresponds to the stretching vibration of O-H bond,free of intra-molecule coupling due to large differences in the atomicmass of H and D74. Two O-H bond length of (1.008 ± 0.002) Å and(1.004 ±0.001) Å given by neutron diffraction at 123 K were thermallycorrected by +0.008Å and averaged to (1.014 ±0.0011) Å75.For ice II, two strong Raman peaks were found at 3194 and3314 cm−1 at 77 K76,77. Neutron diffraction of D2O ice II gives four O-Dbond lengths of (1.014 ± 0.020) Å, (0.975 ± 0.020) Å, (0.956 Å ± 0.020)and (0.937 ±0.020) Å on two symmetrically non-degenerated watermolecules at 110K, with isotropic temperature parameter (B) of 4.5 Å−2for both O and D78. The O-D bond lengths were thermal-correctedusing B, and upshifted further by 3% for the D-H conversion. The finalO-H bond lengths in ice II are (1.102 ± 0.021) Å, (1.064 ±0.021) Å,(1.046 Å ±0.020), and (1.028 ±0.021) Å. The low frequency peak at3194 cm−1 is assigned to the local stretching of the longestO-H bond of1.102 Å. The stretching vibration can be coupled both intra- andintermolecularly among the remaining three O-H bonds due to theirclose bond lengths, resulting in uncertainty in vibration frequencyassignment. Hence, we excluded these three data points and onlyfocused on the longest bond.Ice VI contains two types of O4 octahedra structure inter-penetrating into eachother, revealed by the neutrondiffractionof D2Oice79. For the first O4 octahedra network, all proton positions areequivalent, yielding one O-D bond length (0.986 ± 0.048) Å, indicatinga strong intramolecular coupling in stretching vibration. The secondO4 octahedra structure has three O-D bond lengths at(0.937 ± 0.055) Å, (0.976 ±0.051) Å, and (0.939 ± 0.032) Å with multi-plicity 1, 1, and 2, yielding an average of (0.951 ± 0.027) Å. These twobonds are thermal-corrected to (0.996 ±0.048) Å and(0.961 ± 0.027) Å using B = 2.99 Å−2 for oxygen and 3.80 forhydrogen79, and further upshifted to (1.026 ±0.049) Å and(0.993 ± 0.028) Å for the D-H conversion. These two bonds are asso-ciated with Raman stretching modes at 3390 and 3329 cm−180.For ice VIII, only one bond length is given at (0.969 ±0.007) Å inD2O ice at 10 K81, which is thermal-corrected to (0.981 ± 0.007) Å usingB =0.58Å−2 for O and 1.54 for D, and further upshifted to(1.010 ±0.009) Åbecause ofD-H conversion. TheO-HRamanvibrationfrequency is taken as the average of three intramolecular-couplednormal modes locating at 3477.4 cm−1 (ν1B1g), 3447.6 cm−1 (ν3Eg) and3358.8 cm−1 (ν1A1g) at 100K82, which is 3427.9 cm−1.Ice IX is a proton-ordered phase. Neutron diffraction data on D2Oice shows that protons are preferably arranged in one configuration(referred to as ‘major site’), with a nonzero occupancy on the otherconfiguration (referred to as ‘minor site’). The O-D lengths for majorand minor configurations are (0.976 ±0.003) Å and (1.02 ±0.06) Å.Thermal-corrected length is reported for only the major site at(0.983 ±0.004) Å83. This value is further upshifted to (1.012 ± 0.004) Åfor D-H conversion. Because of the low occupancy, we exclude theminor configuration bond length and only use that of the major con-figuration. Two intense Raman peaks of O-H stretching are found atArticle https://doi.org/10.1038/s41467-025-58608-6Nature Communications | (2025) 16:3447 9www.nature.com/naturecommunications3280.7 cm−1 and 3160.4 cm−1 at 100K84, with the latter assigned to thevibration of O-H bond on major sites.For ice XVII, two O-H stretching Raman peaks at 3106.3 cm−1 and3231.8 cm−1 are reported in ref. 85. These are regarded as decoupledvibration because their separation is >100 cm−1. Two O-D bondlengths are given by neutron diffraction on D2O ice at(1.023 ± 0.007) Å and (1.006 ± 0.004) Å, which are thermal-correctedto (1.044 ± 0.007) Å and (1.018 ± 0.004) Å using B = 2.37 Å−2 for O,4.05 Å−2 for D1 and 3.34 Å−2 for D2, respectively86. These values arefurther upshifted to (1.075 ± 0.007) Å and (1.049 ± 0.004) Å for D-Hconversion.For CaSO4·2H2O (gypsum), the Raman frequency of O-H stretch-ing is taken from this work, which is 3406 cm−1 and 3484 cm−1, mea-sured at room temperature. The O-H bond lengths are taken fromref. 20, (0.961 ± 0.006) Å and (0.948 ± 0.004) Å at 320K, which arethermal-corrected to (0.984 ± 0.006) Å and (0.966 ± 0.004) Å usingB = 1.749Å−2 for O, 3.520Å−2 for D1 and 3.074 Å−2 for D2, respectively.These values are further upshifted to (1.014 ± 0.006) Å and(0.995 ±0.004) Å for D-H conversion.For Li2SO4·H2O, two thermal-corrected O-H bond lengths are(1.004 ±0.003) Å and (0.997 ±0.002) Å87. We associate these twobonds with two Raman peaks of O-H stretching at 3439 cm−1 and3480 cm−1, respectively88.CuSO4·5H2O contains five water molecules per unit cell, amongwhich four are coordinated to Cu2+, and the remaining one is far awayfrom Cu2+ thus non-coordinated. To exclude the interactions withmetal ions, we only focus on the non-coordinated water molecules inthis study. CuSO4·5H2O exhibits multiple water stretching peaks, andthe modes above 3300 cm−1 is assigned to the stretching of non-coordinated water molecules, at 3360 cm−1 and 3477 cm−189,90. The O-Hbond lengths of the non-coordinated water molecules are(0.978 ± 0.005) Å and (0.936 ± 0.011) Å91. These values are correctedfor thermal motion using anisotropic temperature parameters tensor(βij) to (1.005 ±0.005) Å and (0.951 ± 0.011) Å. These two bonds areassociated with 3360 cm−1 and 3477 cm−1 Raman peaks, respectively.For Ba(ClO3)2·H2O, a single thermal-corrected O-H bond length isreported at (0.958± 0.011) Å92, suggesting an intramolecularly coupledvibration. Thus the vibration frequency is taken as the averagebetween ν1 (3512 cm−1) and ν3 (3582 cm−1) vibration modes at 90K93,that is 3547 cm−1.For LiClO4·3H2O, the O-H bond length is taken as the averageamong three refinement series, (0.993 ± 0.004) Å, with uncertaintycoming from the standard deviation of the three series94. The valuesprovided in the reference have been thermal corrected during refine-ment. This bond length is associated with a sharp Raman peak in waterstretching region, 3553 cm−195.For BaCl2·2H2O, four O-H bond lengths are given by neutron dif-fraction measurements at (0.953± 0.003) Å, (0.959±0.003) Å,(0.960 ±0.003) Å, and (0.972 ±0.004) Å96. These values are thermalcorrected to (1.0359± 0.003) Å, (1.0364 ± 0.003) Å, (1.0717 ± 0.003) Åand (1.077 ±0.004) Å using βij. and assigned to four intermolecularlyuncoupled Raman modes at 3456, 3352, 3318, and 3303 cm−197.The assignment and references are given in SupplementaryTable 1, with all the data points plotted in Fig. 3b in the main text. Thedata points above were then fitted using a linear function with both Xand Y errors, yielding a slope of (7.83 ± 1.26) × 10−19C and intercept of(0.957 ±0.001) Å, corresponding to ∂p∂d and d 0ð Þ in Eq. 3 in the maintext. The uncertainty corresponds to 2σ. The reduced χ2 for this fit-ting is 4.96.Dipole moments vs O-H bond lengths in Fig. 3cIn the gas phase, equilibrium dipole moment of H2O molecule is1.855 ± 0.001 D, given by Stark effect measurements29. The H −O −Hangle in gas phase is 104.48° and the O −H bond length is 0.9578 Å, asdetermined by fitting experimental rotational energy levels with therotation-vibration Hamiltonian26. This gives a dipole moment of theO −H bond pOH = 1.514 ± 0.001 D.Effective dipole moment of liquid water at ambient conditions is2.9 ± 0.6 D, estimated from 0.5e (±20%) charge transfer along eachO −H bond, using synchrotron x-ray diffraction measurements32. TheH −O −H angle and the O −H bond length are 103.924° and 0.984Å,obtained from a joint refinement of X-ray and neutron diffractiondata98. No iso- or anisotropic temperature parameters available, so nocorrection was made to the O-H bond length. This gives dipolemoment of O-H bond pOH= 2.4 ± 0.5 D.Ice Ih has a dipole moment of 3.09 ±0.04D, calculated usingmultipole iterative method, with error extracted from Fig. 1 in theref. 33 The O-H bond length and H-O-H angle are 1.014 Å and 109.21°from neutron diffraction data at 60K. The bond length is furtherthermal corrected to 1.0125 Å75. The pOH is calculated to be2.67 ± 0.04D.Polarizability and dielectric behavior of confined waterWe consider a dielectric with εr = ε/ε0 consisting of polarizable dipoleswith the polarizability along the field direction αc,i and the numberdensity of dipoles Ni. The polarization under the electric field E isP = ðεr � 1Þε0E, which equals to the sum of the polarization from alldipoles, P =Piαc, iENi, here i denotes different types of dipoles. Thisleads to Eq. 8 in the main text, the dielectric constant of a confinedwater system.From Eqs. 1, 2 and 6 in the main text, polarizability can be writtenin terms of total electric field E along O-H bond:αðEÞ=∂p∂d� �2k02k0 � ∂2p∂d2 E� �3 ð14ÞTaking gypsum as an example of nano-confined water system,the polarizability of the four types of O-H dipoles in gypsum aredoubly degenerated and dictated by the local EHB field, 5.33 and3.82 V nm−1 for O-HA and O-HB, respectively. In the presence ofexternal electric field, the four dipoles have four different local fields,and thus the structure symmetry is broken, resulting in two polar-izabilities splitting towards two directions (Supplementary Fig. 5c).By adding up the four field-dependent polarizabilities according toEqs. 7 and 8 in the main text, dielectric constant εr = 6. 6 ± 0.7 wasobtained for water in gypsum. The calculated field-tunable εr(E) isplotted in Supplementary Fig. 5d, which presents a modulation of<0.1% due to the field-tuning effect on confined water layers ingypsum.To observe the field-tunable εr(E) in gypsum, we performedcapacitance measurements in one of our devices, using capacitancebridge (Supplementary Fig. 5a). In the measurement, a graphite/gyp-sum/graphite capacitor fabricated on SiO2 (290 nm)/Si substrate wasplaced in vacuum (<10−6mbar). The capacitance of the device wasmeasured using AH 2700 capacitance bridge with an ac excitation of1 V at 1.2 kHz.The equivalent circuit is sketched in Supplementary Fig. 5a.The experimental capacitance measured by the bridge is a combina-tion of four capacitors with three contributions: capacitance of thegypsum (Cg), quantum capacitance of top and bottom graphite elec-trodes (Cq), and parasitic capacitance from the measurement setup(Cp). From capacitors model we can fit the experimental curve using:C =Cp +11=Cq1 + 1=Cq2 + 1=Cgð15ÞField-independent Cp is a fitting parameter with constant value. εr,which is embedded in Cg, is another fitting parameter. 1/Cq1 and 1/Cq2are function of carrier density on graphite electrodes and are given byArticle https://doi.org/10.1038/s41467-025-58608-6Nature Communications | (2025) 16:3447 10www.nature.com/naturecommunicationsref. 99, andweobtainCp ≈ 319 fF and εr ≈ 3.68.WithCp subtracted fromthe experimental curve, the modulation of capacitance is ~0.9%, muchhigher than the 0.1% modulation expected for field-tunable polariz-ability of confined water. Thus, it is concluded that the effect of field-tunable polarizability of water molecules is hidden in the graphitequantum capacitance.Gypsum is a HBH consisting of confined water layers and CaSO4layers, which canbemodeled as two capacitors in series. The total εr ofgypsum has been determined above. The capacitance of confinedwater layers can be evaluated using number density of confined watermolecules and Eq. 8 in the main text. Gypsum (CaSO4·2H2O) hasmolecular weight of 172.17 g·mol−1 and mass density 2.32 g·cm−3100,yielding a number density of 8.11 × 1027m−3 of gypsum, and1.62 × 1028m−3 of H2O. Considering that the water molecules are con-fined between CaSO4 layers, with confining space taking up 40% of thegypsum volume (confining length measured between OS and OS’, theoxygen atoms in SO4 groups fromadjacent layers), the number densityof water molecules in the confining space is corrected to4.06 × 1028 m−3. Equation 8 in the main text produces εr = 6.6 ± 0.7 forconfined 2D water layers in gypsum. Using capacitor-in-series model,εr = 2.8 ± 0.5 can be obtained for CaSO4 molecular sheets. This valuematcheswell the reported εr for anhydrite (CaSO4), which is 2.45 ± 0.52for γ-phase and 2.70 ± 0.02 for β-phase43.Our model finds that water sheets confined in gypsum presents alower εr compared to its bulk form due to the restricted rotationalfreedom. Such behavior is similar to water confined in carbon nano-tubes (CNT). Simulations has shown that in CNT-confined water, axialcomponent of εr increases to above 100 with a relaxation time longerthan bulk water, while the εr of water normal to the carbon wallsdecreases drastically (to ~6), with the dipolar relaxation time severelysuppressed101,102. The small normal component of εr is attributed to therestricted dipole rotation, caused by the strong orientation of watermolecules near walls due to induced surface charge103. A similardecrease in εr was also found at the water—C60 fullerenes interface,which is attributed to the electrostatic interaction between water andsurface charge on C60104,105.Dielectric constant of liquid waterWe consider the dielectric response of liquid water coming from twoparts: dipole reorientation, as well as atomic and electronic polariza-tion. The former is given by Langevin-Debye equation106:PdipoleD E=N0p2E3kTð16ÞwhereN0 and p are the number density and dipolemoment of the H2Omolecule, k is Boltzmann constant, T the temperature and E is theelectric field. The number density of H2O molecule is 3.35 × 1028m−3,using the mass density of liquid water. The latter, atomic and elec-tronic polarization, is given by considering individual O-H dipolemoments of number density 2N0 and dipole moment p, withorientation given by Boltzmann distribution of p in external E-field.The induced dipole moment per O-H bond is αE(cos θ)2, where θ is theangle between p and E (defined in Fig. 1b, main text). Integration overall 3D space gives the polarization density of induced dipole:Pa+ e =2αEN0R�11 y2exydyR�11 exydyð17Þwith y = cos θ, and x = pE/kT. Using Taylor expansion on the inte-gral term and keeping to quadratic term, Eq. 17 reduces toPa+ e =2αEN013+145pEkT� �2" #ð18ÞAdding up the two contributions and using P = ðεr � 1Þε0E, weobtain dielectric constant of liquid water:εr = 1 +N03ε0p2kT+2α +2α15pEkT� �2" #ð19Þwhere p and α are be given by our model from Eqs. 6 and 7 in themain text.Although our model is derived from individual water moleculesconfined in gypsum, it successfully reproduced the dipole moment ofliquidwater, 2.66D, where water dipoles are free to rotate and interactwith each other. Using this result and Eq. 7 in the main text, a staticdielectric constant εr = 27 ± 2 was obtained for liquid water. This issignificantly lower than experimental value εr = 79 due to the complexhydrogen bonding network of water clusters in liquid water, whichenhances the interaction between water molecules and thus alters thecooperativity of water molecules with HBs and affects the dielectricresponse. A conventional approach to include such correlation effectis to multiply p2/kT term in Eq. 19 with a dipole correlation factor Gk44.To get experimental εr = 79, Gk of 3.15 should be taken, in the range ofreported Gk, from 2.72 to 3.70 which depends on the different modelsused for the simulation45. With bond length taken as 0.984 Å98, ourmodel yields a dielectric constant value of 78.7 ± 7.2 at E =0V nm-1 and78.8 ± 7.2 at E =0.5 V nm−1, with the electronic polarization contributesto only 3% of the dipole reorientation.From the calculations above, it is seen that in bulk water, dipolereorientation contributes over 97% to the dielectric response. Previousstudies have shown that in confined liquidwater, particularly when theconfinement is less than 2 nm, the perpendicular dielectric constantdrastically decreases to below 1041,42,49. This significant reduction is aresult of restricted rotational freedom. To quantify the contribution ofelectronic and atomic polarization, which is highly anisotropic anddependent on the water molecule orientation, we calculated thedielectric constant resulting from this effect only. As illustrated inSupplementary Fig. 6, the dielectric constant has a maximum of 13from non-reorientation polarization. Given the experimental εr = 79,dipole reorientation consistently contributes more than 84% of thedielectric response.Consequently, in scenarios where confinement restricts dipolereorientation, and the dielectric response is entirely due to electronicand atomic polarization, water consistently exhibits a low εr.Water-carbon interfaceIn the case of single water molecules trapped inside a C60 fullerene,both the symmetric v1 and antisymmetric v3 stretching modes aresignificantly red-shifted compared to free H2O, to 3573 cm−1 and3659.6 cm−1, respectively. In a harmonic approximation, this corre-sponds to a force constant k ≈ 725Nm−1. Using Eqs. 3 and 4 in themaintext, we obtain a local E-field EHB = 1.7 ± 0.2 V nm−1, O-H bond lengthdOH =0.974 ± 0.004Å, and dipole moment (using Eq. 6 in the maintext) pH2O = 2.37 ± 0.14D. With our dipole-in-E-field model, the calcu-lated interaction strength UHB = 1.5 ± 0.2 kcalmol−1 per bond(66 ± 9meV), which is comparable to typical weak HBs. Similar valuesare obtained for water confined in single-walled carbon nanotubes,with stretching frequency at 3569 and 3640 cm−1, yieldingk ≈ 720Nm−153. However, this picture is likely complicated by theinterplay between water-water intermolecular HBs and water-carbon π HBs.At the water-graphene interface, water molecules exhibit a peakaround 3600 cm−1, which was assigned to the stretching vibration ofthe ‘dangling’ O-H dipole pointing towards the graphene film54. Thereported vibration frequency of νOH = 3616 cm−1 for the dangling O-Hdipole, using our approach, corresponds to k = 725.1Nm−1,EHB = 1.7 ± 0.2 V nm−1, dOH = 0.975 ± 0.004Å, pH2O = 2.37 ± 0.14D, andUHB = 1.54 ± 0.20 kcalmol−1 per bond (67 ± 9meV). Our model predictsArticle https://doi.org/10.1038/s41467-025-58608-6Nature Communications | (2025) 16:3447 11www.nature.com/naturecommunicationsthe electric field exerted by theπ electrons on the water O-H dipoles isaround 1.7 V nm−1. Using only the stretching frequency of confinedwater, our results directly reproduce key properties of confinedwater,including HB strength, local E-field, O-H bond length, and dipolemoment.This peak can also be tuned by applying an external electric fieldin the electricdouble layer regionneargraphene surface. In the vicinityof graphene, water O-H dipoles experience an E-field exerted from πelectrons of carbon atoms. This interaction can also be considered aπ-hydrogen bond (π-HB),which redshifts the stretching frequency of theO-H dipole compared to that of a free water molecule. In the presenceof Eext, the interaction is tuned, which can be seen in the shift of thestretching frequency of dangling O-H bond on the graphene surface.Reference 54 reported a sum frequency vibrational spectroscopy(SFVS) study on the water-graphene interface. As an independentvalidation of ourmodel, we used the data points from Fig. 4b in ref. 54.The vibration frequency of danglingO-H peak on y-axis is converted toforce constant k using harmonic approximation. The interfacial elec-tric field Es on x-axis corresponds to the Eext in our study. Etot for eachpoint is obtained using Eq. 5 inmain text, with EHB given by Eq. 4 (maintext) using k at Eext = 0, taken fromFig. 4b in ref. 54. Since the danglingO-H dipole points towards graphene film, cos ϑ in main text Eq. (5) istaken as 1. The obtained results match perfectly into our model, asplotted in Supplementary Fig. 7.Dielectric constants of water ice systemsDue to the presence of defects in proton-disordered ices which causedipole reorientation and, thus, an anomalously large εr, we chose theproton-ordered ices to study the static εr. Hence, only atomic andelectronic polarization contributions are considered. Ice XI is proton-ordered form of ice Ih, with εr = 2.9 taken from Fig. 6 in ref. 107. Theauthors cooled down the ordinary ice (ice Ih) and measured εr at 77 K,where the ice Ih–XI phase transition is expected to occur. Theobserved small and non-dispersive εr indicates a proton-ordered formof ordinary ice. The overall static εr of proton-ordered ice II is around4.2, measured from Fig. 6 in ref. 47, although a small increase wasobserved with frequency (not evaluated). Ice IX is proton-orderedformof ice III. The εr of ice IX is taken as4 according to Fig. 3 in ref. 108,where static εr presented a sharp drop during the ice III–IX phasetransition.We examined the static dielectric constant along [001], [111] and[010] directions for these three phases of ices, respectively, axesdirection defined in refs. 109–111. The three confined water systemspossess 2, 4, and 6 types of non-equivalent O-H dipoles with respect tothe assessing axis, with polarizability marked in Fig. 3d in main text.Using the neutron diffraction data of O-H bond lengths reported inrefs. 75,78,112 (the bond length of ice III is thermal corrected to 1.005 Åfollowing the approach in Methods, ‘Thermal corrections for neutronstructural data’) and density reported in refs. 110,113, the dielectricconstants calculated using Eq. 8 in themain text, are 3.8 ± 0.5, 5.9 ± 0.8,and 4.7 ± 0.4 for structure of ice XI, II and IX, respectively. Our resultsshow reasonable agreement with the reported εr for ice XI, II, and IX,which are 2.9107, 4.247 and 4108.All the estimated εr for water systems are summarized in Sup-plementary Table 2.Data availabilityThe SourceDatafiles for all data presented ingraphswithin the Figuresused in this study are available in the Zenododatabase: https://doi.org/10.5281/zenodo.15019875.References1. Gordon, M. S. & Jensen, J. H. Understanding the hydrogen bondusing quantum chemistry. Acc. Chem. Res. 29, 536–543 (1996).2. Grabowski, S. J. Hydrogen Bonding - New Insights (Springer Dor-drecht, 2006).3. Li, X.-Z., Walker, B. & Michaelides, A. Quantum nature of thehydrogen bond. Proc. Natl Acad. Sci. USA 108, 6369–6373 (2011).4. Herschlag, D. & Pinney, M. M. Hydrogen bonds: simple after all?Biochemistry 57, 3338–3352 (2018).5. Arunan, E. et al. Definition of the hydrogen bond (IUPAC Recom-mendations 2011). Pure Appl. Chem. 83, 1637–1641 (2011).6. Arunan, E. et al. Defining the hydrogen bond: an account (IUPACTechnical Report). Pure Appl. Chem. 83, 1619–1636 (2011).7. Dereka, B. et al. Crossover from hydrogen to chemical bonding.Science 371, 160–164 (2021).8. Ceriotti, M., Cuny, J., Parrinello, M. & Manolopoulos, D. E. Nuclearquantum effects and hydrogen bond fluctuations in water. Proc.Natl. Acad. Sci. USA 110, 15591–15596 (2013).9. Van Hoozen, B. L. Jr. & Petersen, P. B. Vibrational tug-of-war: thepK(A) dependence of the broad vibrational features of stronglyhydrogen-bonded carboxylic acids. J. Chem. Phys. 148, 134309(2018).10. Cleland, W. W. & Kreevoy, M. M. Low-barrier hydrogen bonds andenzymic catalysis. Science 264, 1887–1890 (1994).11. Steiner, T. The hydrogenbond in the solid state.Angew.Chem. Int.Ed. Engl. 41, 49–76 (2002).12. Carrasco, J., Hodgson, A. & Michaelides, A. A molecular per-spective of water at metal interfaces. Nat. Mater. 11, 667–674(2012).13. Pauling, L. The shared-electron chemical bond. Proc. Natl. Acad.Sci. USA 14, 359–362 (1928).14. Yamabe, S. &Morokuma, K. Molecular orbital studies of hydrogenbonds. IX. Electron distribution analysis. J. Am. Chem. Soc. 97,4458–4465 (1975).15. Paoloni, L. Nature of the Hydrogen Bond. J. Chem. Phys. 30,1045–1058 (1959).16. Kolesov, B. A. Hydrogen bonds: Raman spectroscopic study. Int. J.Mol. Sci. 22, 5380 (2021).17. Hush, N. S. &Williams, M. L. Carbonmonoxide bond length, forceconstant and infrared intensity variations in strong electric fields:valence-shell calculations, with applications to properties ofadsorbed and complexed CO. J. Mol. Spectrosc. 50, 349–368(1974).18. Cole, W. F. & Lancucki, C. J. Symmetry of SO4 Ion in gypsum.Nat.Phys. Sci. 238, 95–96 (1972).19. Cole, W. F. & Lancucki, C. J. Hydrogen bonding in gypsum. Nat.Phys. Sci. 242, 104–105 (1973).20. Schofield, P. F., Knight, K. S. & Stretton, I. C. Thermal expansion ofgypsum investigated by neutron powder diffraction. Am. Miner-alo. 81, 847–851 (1996).21. Ohno, K., Okimura, M., Akai, N. & Katsumoto, Y. The effect ofcooperative hydrogen bonding on the OH stretching-band shiftfor water clusters studied by matrix-isolation infrared spectro-scopy and density functional theory. Phys. Chem. Chem. Phys. 7,3005–3014 (2005).22. Frisenda, R. et al. Recent progress in the assembly of nanodevicesand vanderWaals heterostructures bydeterministic placement of2D materials. Chem. Soc. Rev. 47, 53–68 (2018).23. Couty, R., Velde, B. & Besson, J. M. Raman spectra of gypsumunder pressure. Phys. Chem. Miner. 10, 89–93 (1983).24. Seki, T. et al. The bending mode of water: a powerful probe forhydrogen bond structure of aqueous systems. J. Phys. Chem. Lett.11, 8459–8469 (2020).25. Furtenbacher, T., Tóbiás, R., Tennyson, J., Polyansky, O. L. &Császár, A. G. W2020: a database of validated rovibrationalexperimental transitions and empirical energy levels of H216O. J.Phys. Chem. Ref. Data 49, 033101 (2020).Article https://doi.org/10.1038/s41467-025-58608-6Nature Communications | (2025) 16:3447 12https://doi.org/10.5281/zenodo.15019875https://doi.org/10.5281/zenodo.15019875www.nature.com/naturecommunications26. Hoy, A. R. & Bunker, P. R. A precise solution of the rotation bendingSchrödinger equation for a triatomic molecule with application tothe water molecule. J. Mol. Spectrosc. 74, 1–8 (1979).27. Csaszar, A. G. et al. On equilibrium structures of the water mole-cule. J. Chem. Phys. 122, 214305 (2005).28. Beers, Y. & Klein, G. P. The stark splitting of millimeter wavetransitions of water. J. Res. Natl. Bur. Stand A Phys. Chem. 76A,521–528 (1972).29. Clough, S. A., Beers, Y., Klein, G. P. & Rothman, L. S. Dipolemoment of water from Stark measurements of H2O, HDO, andD2O. J. Chem. Phys. 59, 2254–2259 (1973).30. Lovas, F. J. Microwave spectral tables II. Triatomic molecules. J.Phys. Chem. Ref. Data 7, 1445–1750 (1978).31. Lodi, L. et al. A new ab initio ground-state dipole moment surfacefor the water molecule. J. Chem. Phys. 128, 044304 (2008).32. Badyal, Y. et al. Electron distribution in water. J. Chem. Phys. 112,9206–9208 (2000).33. Batista, E. R., Xantheas, S. S. & Jónsson, H. Molecular multipolemoments of water molecules in ice Ih. J. Chem. Phys. 109,4546–4551 (1998).34. Bjerrum, N. Structure and properties of ice. Science 115, 385–390(1952).35. Ishii, F., Terada, K. &Miura, S. First-principles studyof spontaneouspolarisation and water dipole moment in ferroelectric ice XI.Mol.Simul. 38, 369–372 (2012).36. Carey, D. M. & Korenowski, G. M. Measurement of the Ramanspectrum of liquid water. J. Chem. Phys. 108, 2669–2675 (1998).37. Walrafen, G. Raman spectral studies of the effects of temperatureon water structure. J. Chem. Phys. 47, 114–126 (1967).38. Lide, D. R. CRC Handbook of Chemistry and Physics (CRCPress, 2004).39. Whalley, E. Energies of the phases of ice at zero temperature andpressure. J. Chem. Phys. 81, 4087–4092 (1984).40. Nanayakkara, S., Tao, Y. & Kraka, E. Capturing individual hydrogenbond strengths in ices via periodic local vibrational mode theory:beyond the lattice energy picture. J. Chem. Theory Comput. 18,562–579 (2021).41. Schlaich, A., Knapp, E. W. & Netz, R. R. Water dielectric effects inplanar confinement. Phys. Rev. Lett. 117, 048001 (2016).42. Loche, P., Ayaz, C., Wolde-Kidan, A., Schlaich, A. & Netz, R. R.Universal and nonuniversal aspects of electrostatics in aqueousnanoconfinement. J. Phys. Chem. B 124, 4365–4371 (2020).43. López-Buendía, A. M., García-Baños, B., Urquiola, M. M., Catalá-Civera, J. M. & Penaranda-Foix, F. L. Evidence of a new phase ingypsum–anhydrite transformations under microwave heating byin situ dielectric analysis and Raman spectroscopy. Phys. Chem.Chem. Phys. 22, 27713–27723 (2020).44. Caillol, J. Asymptotic behavior of the pair‐correlation function of apolar liquid. J. Chem. Phys. 96, 7039–7053 (1992).45. Yu, H., Hansson, T. & van Gunsteren, W. F. Development of asimple, self-consistent polarizable model for liquid water. J.Chem. Phys. 118, 221–234 (2003).46. Gränicher, H., Jaccard, C., Scherrer, P. & Steinemann, A. Dielectricrelaxation and the electrical conductivity of ice crystals. Discuss.Faraday Soc. 23, 50–62 (1957).47. Wilson, G., Chan, R., Davidson, D. & Whalley, E. Dielectric prop-erties of ices II, III, V, and VI. J. Chem. Phys. 43, 2384–2391 (1965).48. Bramwell, S. T. Ferroelectric ice. Nature 397, 212–213 (1999).49. Fumagalli, L. et al. Anomalously low dielectric constant of con-fined water. Science 360, 1339–1342 (2018).50. Michaelides, A. et al. Origin of dielectric polarization suppressionin confined water from first principles. Chem. Sci. 15, 516–527(2023).51. Bocquet, L. Nanofluidics coming of age. Nat. Mater. 19, 254–256(2020).52. Shugai, A. et al. Infrared spectroscopy of an endohedral water infullerene. J. Chem. Phys. 154, 124311 (2021).53. Dalla Bernardina, S. et al. Water in carbon nanotubes: the peculiarhydrogenbondnetwork revealed by infrared spectroscopy. J. Am.Chem. Soc. 138, 10437–10443 (2016).54. Xu, Y.,Ma, Y.-B.,Gu, F., Yang,S.-S. &Tian,C.-S. Structure evolutionat the gate-tunable suspended graphene–water interface. Nature621, 506–510 (2023).55. Keutsch, F. N. & Saykally, R. J. Water clusters: untangling themysteries of the liquid, one molecule at a time. Proc. Natl. Acad.Sci. USA 98, 10533–10540 (2001).56. Zhang, X. et al. Large-scale 2D heterostructures from hydrogen-bonded organic frameworks and graphene with distinct Dirac andflat bands. Nat. Commun. 15, 5934 (2024).57. Daldrop, J. O. et al. Orientation of non-spherical protonated waterclusters revealed by infrared absorption dichroism.Nat. Commun.9, 311 (2018).58. Martin, L. W. & Rappe, A. M. Thin-film ferroelectric materials andtheir applications. Nat. Rev. Mater. 2, 1–14 (2016).59. Zhang, Z., Ye, Y., Xiang, S. & Chen, B. Exploring multifunctionalhydrogen-bonded organic framework materials. Acc. Chem. Res.55, 3752–3766 (2022).60. Shi, B. et al. Short hydrogen-bond network confined on COF sur-faces enables ultrahigh proton conductivity. Nat. Commun. 13,6666 (2022).61. Lin, Z. J., Mahammed, S. A. R., Liu, T. F. & Cao, R. Multifunctionalporous hydrogen-bonded organic frameworks: current status andfuture perspectives. ACS Cent. Sci. 8, 1589–1608 (2022).62. Pal, S. C., Mukherjee, D., Sahoo, R., Mondal, S. & Das, M. C. Proton-conducting hydrogen-bonded organic frameworks. ACS EnergyLett. 6, 4431–4453 (2021).63. Li, P., Ryder, M. R. & Stoddart, J. F. Hydrogen-bonded organicframeworks: a rising class of porous molecular materials. Acc.Mater. Res. 1, 77–87 (2020).64. Castellanos-Gomez, A. et al. Deterministic transfer of two-dimensional materials by all-dry viscoelastic stamping. 2D Mater.1, 011002 (2014).65. Pizzocchero, F. et al. The hot pick-up technique for batch assem-bly of van der Waals heterostructures. Nat. Commun. 7, 11894(2016).66. Chio,C. H., Sharma, S. K. &Muenow,D.W.Micro-Raman studies ofgypsum in the temperature range between 9 K and 373 K. Am.Mineral. 89, 390–395 (2004).67. Flaud, J. M. & Camy-Peyret, C. The interacting states (020), (100),and (001) of H216O. J. Mol. Spectrosc. 51, 142–150 (1974).68. Flaud, J. M. & Camy-Peyret, C. Vibration-rotation intensities inH2O-type molecules application to the 2ν2, ν1, and ν3 bands ofH216O. J. Mol. Spectrosc. 55, 278–310 (1975).69. Murphy, W. F. The rovibrational Raman spectrum of water vapourv1 and v3. Mol. Phys. 36, 727–732 (1978).70. Rahn, L. A. & Greenhalgh, D. A. High-resolution inverse Ramanspectroscopy of the ν1 band ofwater vapor. J. Mol. Spectrosc. 119,11–21 (1986).71. Tennyson, J. et al. IUPAC critical evaluation of therotational–vibrational spectra of water vapor. Part II. J. Quant.Spectrosc. Radiat. Transf. 111, 2160–2184 (2010).72. Busing, W. R. & Levy, H. A. The effect of thermal motion on theestimation of bond lengths from diffraction measurements. ActaCrystallogr. 17, 142–146 (1964).73. Soper, A. K. & Benmore, C. J. Quantumdifferences between heavyand light water. Phys. Rev. Lett. 101, 065502 (2008).74. Scherer, J. R. & Snyder, R. G. Raman intensities of single crystal iceI h. J. Chem. Phys. 67, 4794–4811 (1977).75. Kuhs,W. F. & Lehmann,M. S. The structure of the ice Ih by neutrondiffraction. J. Phys. Chem. 87, 4312–4313 (1983).Article https://doi.org/10.1038/s41467-025-58608-6Nature Communications | (2025) 16:3447 13www.nature.com/naturecommunications76. Taylor, M. &Whalley, E. Raman spectra of ices Ih, Ic, II, III, and V. J.Chem. Phys. 40, 1660–1664 (1964).77. Minceva-Sukarova, B., Sherman, W. & Wilkinson, G. The Ramanspectra of ice (Ih, II, III, V, VI and IX) as functions of pressure andtemperature. J. Phys. C Solid State Phys. 17, 5833 (1984).78. Kamb, B., Hamilton, W. C., LaPlaca, S. J. & Prakash, A. Orderedproton configuration in ice II, from single‐crystal neutron diffrac-tion. J. Chem. Phys. 55, 1934–1945 (1971).79. Kuhs, W. F., Finney, J. L., Vettier, C. & Bliss, D. V. Structure andhydrogen ordering in ices VI, VII, and VIII by neutron powder dif-fraction. J. Chem. Phys. 81, 3612–3623 (1984).80. Thoeny, A. V., Parrichini, I. S., Gasser, T. M. & Loerting, T. Ramanspectroscopy study of the slow order-order transformation ofdeuterium atoms: ice XIX decay and ice XV formation. J. Chem.Phys. 156, 154507 (2022).81. Jorgensen, J. D., Beyerlein, R., Watanabe, N. & Worlton, T. Struc-ture of D2O ice VIII from insitu powder neutron diffraction. J.Chem. Phys. 81, 3211–3214 (1984).82. Wong, P. & Whalley, E. Raman spectrum of ice VIII. J. Chem. Phys.64, 2359–2366 (1976).83. La Placa, S. J., Hamilton, W. C., Kamb, B. & Prakash, A. On a nearlyproton‐ordered structure for ice IX. J. Chem. Phys. 58, 567–580(1973).84. Bertie, J. E. & Francis, B. F. Raman spectra of the O–H and O–Dstretching vibrations of ices II and IX to 25 K at atmosphericpressure. J. Chem. Phys. 72, 2213–2221 (1980).85. Del Rosso, L. et al. Cubic ice Ic without stacking defects obtainedfrom ice XVII. Nat. Mater. 19, 663–668 (2020).86. del Rosso, L. et al. Refined structure of metastable ice XVII fromneutron diffraction measurements. J. Phys. Chem. C 120,26955–26959 (2016).87. Lundgren, J. O., Kvick, Å., Karppinen, M., Liminga, R. & Abrahams,S. Neutron diffraction structural study of pyroelectric Li2SO4⋅H2O at 293, 80, and 20. K. J. Chem. Phys. 80, 423–430 (1984).88. Canterford, R. & Ninio, F. Raman study of some water modes inlithium sulphate monohydrate. Solid State Commun. 15,1451–1452 (1974).89. Liu, D. & Ullman, F. Raman spectrum of CuSO4· 5H2O singlecrystal. J. Raman Spectrosc. 22, 525–528 (1991).90. Berger, J. Infrared and Raman spectra of CuSO4, 5H2O; CuSO4,5D2O; and CuSeO4, 5H2O. J. Raman Spectrosc. 5, 103–114 (1976).91. Bacon, G. & Titterton, D. Neutron-diffraction studies of CuSO4·5H2O and CuSO4· 5D2O. Z. F. Kristallographie-Crystalline Mater.141, 330–341 (1975).92. Sikka, S., Momin, S., Rajagopal, H. & Chidambaram, R. Neu-tron‐diffraction refinement of the crystal structure of bariumchlorate monohydrate Ba (ClO3) 2· H2O. J. Chem. Phys. 48,1883–1889 (1968).93. Bertie, J. E., Heyns, A. M. & Oehler, O. The infrared and Ramanspectra of powdered and single-crystal barium chlorate mono-hydrate and-deuterate. Can. J. Chem. 51, 2275–2289 (1973).94. Lundgren, J.-O., Liminga, R. & Tellgren, R. Neutron diffractionrefinement of pyroelectric lithium perchlorate trihydrate. ActaCrystallogr. Sect. B Struct. Crystallogr. Cryst. Chem. 38,15–20 (1982).95. Zhang, Y.-H. & Chan, C. K. Observations of water monomers insupersaturated NaClO4, LiClO4, and Mg (ClO4) 2 droplets usingRaman spectroscopy. J. Phys. Chem. A 107, 5956–5962 (2003).96. Padmanabhan, V., Busing, W. & Levy, H. A. Barium chloride dihy-drate by neutron diffraction. Acta Crystallogr. Sect. B Struct.Crystallogr. Cryst. Chem. 34, 2290–2292 (1978).97. Lutz, H., Pobitschka, W., Frischemeier, B. & Becker, R.-A. LatticeVibration Spectra. XX. Infrared and Raman Spectra of BaCl 2• 2H 2O and BaCl 2• 2D 2. O. Appl. Spectrosc. 32, 541–547 (1978).98. Soper, A. Joint structure refinement of x-ray and neutron diffrac-tion data on disordered materials: application to liquid water. J.Phys. Condens. Matter 19, 335206 (2007).99. Mullan, C. et al. Mixing of moiré-surface and bulk states in gra-phite. Nature 620, 756–761 (2023).100. Robie, R. A. & Bethke, P. M. Molar volumes and densities ofminerals. Unites States Department of the Interior Geological Sur-vey, 4–21 (UNT Digital Library, 1962).101. Lin, Y., Shiomi, J., Maruyama, S. & Amberg, G. Dielectric relaxationof water inside a single-walled carbon nanotube. Phys. Rev. B 80,045419 (2009).102. Loche, P., Ayaz, C., Schlaich, A., Uematsu, Y. & Netz, R. R. Giantaxial dielectric response in water-filled nanotubes and effectiveelectrostatic ion–ion interactions from a tensorial dielectricmodel. J. Phys. Chem. B 123, 10850–10857 (2019).103. Renou, R., Szymczyk,A. &Ghoufi, A. Tunabledielectric constant ofwater at the nanoscale. Phys. Rev. E 91, 032411 (2015).104. Sarhangi, S. M., Waskasi, M. M., Hashemianzadeh, S. M. &Matyushov, D. V. Effective dielectric constant of water at theinterface with charged C60 fullerenes. J. Phys. Chem. B 123,3135–3143 (2019).105. Meier, B. et al. Electrical detection of ortho–para conversion infullerene-encapsulated water. Nat. Commun. 6, 8112 (2015).106. Debye, P. Einige resultate einer kinetischen theorie der isolatoren.Phys. Z. 13, 97 (1912).107. Shimizu, N. & Muramoto, Y. Natural materials adopted for elec-trical and electronic materials. IEEJ Trans. Fundamentals Mater.127, 429–432 (2007).108. Whalley, E., Heath, J. & Davidson, D. Ice IX: an antiferroelectricphase related to ice III. J. Chem. Phys. 48, 2362–2370 (1968).109. Bernal, J. D. & Fowler, R. H. A theory of water and ionic solution,with particular reference to hydrogen and hydroxyl ions. J. Chem.Phys. 1, 515–548 (1933).110. Kamb, B. Ice. ii. a proton-ordered form of ice. Acta Crystallogr. 17,1437–1449 (1964).111. Kamb, B. & Prakash, A. Structure of ice III. Acta Crystallogr. Sect. BStruct. Crystallogr. Cryst. Chem. 24, 1317–1327 (1968).112. Londono, J. D., Kuhs, W. F. & Finney, J. L. Neutron diffraction stu-dies of ices III and IX on under‐pressure and recovered samples. J.Chem. Phys. 98, 4878–4888 (1993).113. Bridgman, P. W. Water, in the liquid and five solid forms, underpressure (American Academy of Arts and Sciences, 1913).AcknowledgementsThis research was supported by the European Research Council (ERC)under the European Union’s Horizon 2020 research and innovationprogram (Grant Agreement No. 865590) and the Research Council UK[BB/X003736/1]. Q.Y. acknowledges the funding from Royal SocietyUniversity Research Fellowship URF\R1\221096 and UK Research andInnovation Grant [EP/X017575/1].Author contributionsA.M. and Q.Y. initiated and supervised the project. Z.W. fabricateddevices, performed Raman measurements and analyzed the data withthe help from I.T. and A.N.G. A.B., M.Y., and F.P. performed DFT cal-culations. V.K. and A.N.G. helped with FTIR measurements. P.D.Nperformed scanning near-field microscopy measurements. C.M. car-ried out the low temperature capacitancemeasurements. T.T. and K.W.provided hBN crystals. A.M., Z.W., Q.Y., and K.S.N. contributed towriting the paper. All authors discussed the results and commented onthe paper.Competing interestsThe authors declare no competing interests.Article https://doi.org/10.1038/s41467-025-58608-6Nature Communications | (2025) 16:3447 14www.nature.com/naturecommunicationsAdditional informationSupplementary information The online version containssupplementary material available athttps://doi.org/10.1038/s41467-025-58608-6.Correspondence and requests for materials should be addressed toZiwei Wang, Qian Yang or Artem Mishchenko.Peer review information Nature Communications thanks the anon-ymous reviewer(s) for their contribution to thepeer reviewof thiswork. Apeer review file is available.Reprints and permissions information is available athttp://www.nature.com/reprintsPublisher’s note Springer Nature remains neutral with regard to jur-isdictional claims in published maps and institutional affiliations.Open Access This article is licensed under a Creative CommonsAttribution 4.0 International License, which permits use, sharing,adaptation, distribution and reproduction in any medium or format, aslong as you give appropriate credit to the original author(s) and thesource, provide a link to the Creative Commons licence, and indicate ifchanges were made. The images or other third party material in thisarticle are included in the article's Creative Commons licence, unlessindicated otherwise in a credit line to the material. If material is notincluded in the article's Creative Commons licence and your intendeduse is not permitted by statutory regulation or exceeds the permitteduse, you will need to obtain permission directly from the copyrightholder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.© The Author(s) 2025Article https://doi.org/10.1038/s41467-025-58608-6Nature Communications | (2025) 16:3447 15https://doi.org/10.1038/s41467-025-58608-6http://www.nature.com/reprintshttp://creativecommons.org/licenses/by/4.0/http://creativecommons.org/licenses/by/4.0/www.nature.com/naturecommunications Quantifying hydrogen bonding using electrically tunable nanoconfined water Results Dipole-in-E-field model of hydrogen bonding Experimental realization of the model in gypsum Electrically tunable stretching vibration modes Fitting dipole-in-E-field model with experimental data Dipole-in-E-field model for hydrogen bond energy Predicting the dielectric properties of HB systems Beyond conventional HB systems Discussion Methods Full form of the dipole-in-E-field approximation Device fabrication Vibrational spectroscopy measurements Data analysis of Raman spectra Decoupling of the stretching frequencies of O-H bonds in crystalline water of gypsum Infrared spectroscopy of water O-HA stretching modes in gypsum Water bending mode in gypsum Fitting of k vs E relation in Fig. 3a and local fields at dipole O-HA and O-HB Thermal corrections for neutron structural data Data analysis for bond length vs E/k relation in Fig. 3b Dipole moments vs O-H bond lengths in Fig. 3c Polarizability and dielectric behavior of confined water Dielectric constant of liquid water Water-carbon interface Dielectric constants of water ice systems Data availability References Acknowledgements Author contributions Competing interests Additional information