Introduction

For decades intense research activities tackle the question whether water molecules with their rather strong dipole moment of p0 = 1.85 Debye can condense into a ferroelectrically or antiferroelectrically ordered state. In pure liquid water or in H2O ice no such ordering occurs under ambient conditions or during experimentally available time scales because a complex time-dependent tetrahedral molecular hydrogen-bonded network emerges governed by Pauling’s ice rules1,2,3. Proton ordered phases of water ice can be induced by defects introduced via doping with impurities or via irradiation3,4. It is suggested that the so-called water ferroelectricity can exist in various natural systems where the role of H-bonds for intermolecular coupling is diminished due to a preferred orientation of water molecules imposed by the host frameworks, thus creating conditions favorable for mutual aligning the H2O molecular dipoles. In practice, this can be realized when water molecules are adsorbed by two-dimensional interfaces or confined in nanosized channels or cages. Such kind of ferroelectricity of confined water molecules is discussed to be of central importance for a wide range of phenomena in geology, mineralogy, meteorology, soil chemistry, etc. Of special interest is the potential ordering of water molecules in biological systems, where the H2O molecules are found in fully or partially confined states in cells, membrane channels, and proteins hydration shells5,6,7,8, which can influence the functioning of living organisms substantially.

Various types of dipolar orderings and other exciting properties are predicted based on theoretical analyses or computer simulations of water molecules in nanometer-thick layers on substrates or at the interface9,10,11,12,13, between two graphene layers14,15,16,17, for H2O chains within carbon nanotubes18,19,20,21,22,23, for networks of H2O molecules in nanovoids, e.g., in fullerenes24,25,26 or within protein hydration shells27. Besides fundamental aspects, the interest in such systems is fueled by the perspective that ordered water dipoles can find practical applications in nanoelectronic devices28,29,30,31,32,33. It turned out, however, that it is not so trivial to experimentally verify the predictions made by theory and modeling. Even to realize the dipolar ordering under laboratory conditions appears to be a challenging task. Corresponding reports in the literature are very rare9,10,11,34,35, and either observe only a tiny fraction of ordered dipoles (≤1%)10, or raise questions and discussions concerning the reliability of the results36,37,38,39,40,41.

An ideal workbench for studying collective effects in ensembles of dipole–dipole coupled water molecules is provided by hydrated dielectric crystals, such as the beryl family. Compared to zeolites and microporous silicates42,43, the special appeal of these crystals is that they contain just separate H2O molecules that are isolated within nanosized voids formed by the lattice ions. Being only weakly linked to the ions and separated by 5–10 Å, the water molecules do not experience hydrogen bonding (interaction length 1–2 Å); nevertheless, they interact via long-range dipole–dipole coupling (interaction length 10–100 Å). This kind of network of interacting water molecules is of broad interest and fundamental importance providing the opportunity to study electric-dipolar systems whose properties are expected to be qualitatively different from those occurring in systems with magnetic moments44,45. In addition, introducing dopants like Na and K ions at bottlenecks that separate nanovoids can lead to a polarization of H2O molecules and the formation of statically polarized so-called water-II molecules. Having also the possibility to encapsulate heavy water molecules (D2O, DHO), one gets a unique “laboratory” for studies of various excitations, phases and phase transitions in the water dipolar network. Very recently it was reported46,47,48,49,50,51 that H2O molecules in beryl exhibit the tendency towards a macroscopic alignment of their dipoles: a ferroelectric soft mode was observed with the temperature evolution of dielectric strength and frequency following the Curie–Weiss and Cochran laws, respectively. However, due to quantum tunneling of the dipoles in the symmetric six-well localizing potential of the hexagonal beryl lattice52 no phase transition into a macroscopically ordered phase occurs.

In this communication, we present our studies of the dielectric response, pyrocurrent and specific heat of an array of separate H2O molecules confined within ionic matrix of the orthorhombic cordierite crystal lattice, and supplement our experimental data with Density Functional Theory Molecular Dynamics and Monte Carlo simulations. We discover that at T0 ≈ 3 K the water dipolar network undergoes an order–disorder ferroelectric phase transition. In the ground state the water molecules form ferroelectric domains in the ab-plane that order antiferroelectrically along the channel direction.

Results

Material

In the present study, we explore the possibility to order water molecules macroscopically by dipole–dipole interaction. To that end, we consider hydrated crystals with a strongly asymmetric localizing potential; i.e. the presence of a certain crystallographic direction along which the molecular dipoles align and form an ordered phase. Crystalline cordierite (Mg,Fe)2Al4Si5O18 is ideally suited for these requirements as it is structurally similar to beryl, except for the lack of rotational symmetry. Beryl contains Si6O18 rings, while in cordierite two of six silicon atoms are replaced by aluminum, leading to (Si,Al)6O18 rings that are stacked along the c-axis and form channels with 2.5 Å diameter ionic “bottlenecks”, see Fig. 1. As a consequence, the cages in cordierite are anisotropic (5.4 Å along the b-direction and 6.0 Å along the a-axis) and the lattice exhibits orthorhombic symmetry space group Cccm53. There are indications that the H2O dipole moments in cordierite are preferably oriented parallel to the b-axis54. In the present study, slices of well characterized natural hydrous cordierite single crystal were cut along different crystallographic axes in order to measure the dielectric response in all three principal polarizations. We have investigated the electrodynamic properties in a broad frequency range, ν = 1 Hz–3 THz, covering temperatures down to 0.3 K, supplemented by pyroelectric and heat capacity measurements. From the comparison with reference measurements on dehydrated samples, the characteristics of the water molecular subsystem can be extracted.

Fig. 1: Cordierite crystal structure with water molecules within ionic nanopores.
figure 1

a Unit cell of cordierite crystal from X-ray analysis with thermal ellipsoids at 85 K. b The water molecules within cordierite nanopores form two-dimensional triangular lattice within the ab-plane. c One-dimensional chain of confined water molecules along the c-axis.

Dielectric spectra

For the E || a polarization, we find a broad relaxation-like excitation whose peak frequency decreases when the temperature is lowered down to T0 ≈ 3 K and increases on further cooling. Around the same temperature, pronounced anomalies are observed in the dielectric strength of the excitation and of the H2O contribution to the specific heat. We assign this behavior to an order–disorder phase transition within the water molecular network. Density functional theory molecular dynamics (DFT-MD) and Monte Carlo simulations of water molecules indicate that the low-temperature ordered phase contains domains lying in the ab-planes and composed of unidirectional H2O dipoles polarized preferably along the b-axis but ordered antiferroelectrically along the c-axis.

Figure 2 presents the temperature-dependent radiofrequency (RF) spectra of the real ε′(ν) and imaginary ε″(ν) parts of the dielectric permittivity of water molecules in cordierite measured for the polarization E || a. Since the MHz–GHz reflectometric technique does not provide absolute values, the corresponding loss data are separately displayed in Supplementary Fig. 1. For the polarization E || b, the dielectric response is mostly governed by the higher frequency THz modes and does not reveal features that could be indicative of dipolar ordering. For that reason, in the following we focus on the E || a polarization, where strong variations in the spectra are observed on lowering the temperature. We find a pronounced maximum in ε′(T) with a slope above ≈40 K that follows the Curie–Weiss behavior:

$$\varepsilon^{\prime}\left( T \right) = \varepsilon _\infty + C(T - T_{\mathrm{c}})^{ - 1}.$$
(1)
Fig. 2: Temperature evolution of relaxational absorption band.
figure 2

a Temperature dependences of the real part of the dielectric permittivity ε′ of a hydrous cordierite crystal for E || a polarization at different frequencies as indicated. In addition, ε′(T) data for a water-free cordierite crystal (WF, 1.2 kHz) are included. The blue solid line shows the result of a least-square fit with the Curie-Weiss expression, Eq. (1), with C = 390 K and TC = −15 K. b, c Hz–MHz and THz spectra of the real b and imaginary c parts of the dielectric permittivity of the hydrous cordierite crystal measured at different temperatures as indicated. The short-dashed lines represent least-square fits with Eq. (2) of measured spectra shown by full dots. WF denotes the spectra obtained for a water-free crystal (squares). The dotted lines connecting the radiofrequency and terahertz spectra are guides to the eye.

Here, ε is the high-frequency dielectric constant, C is the Curie constant and TC is the Curie temperature. The low-temperature wing of the maximum reveals a pronounced frequency dispersion that is caused by a relaxational band (Fig. 2b, c) signified by a typical relaxation step in ε′(ν) together with a broad peak in ε″(ν). The fact that this feature is completely absent in water-free crystals unambiguously links the relaxation process to the water subsystem in cordierite. Upon cooling down to T0 ≈ 3 K, the band shifts to lower frequencies, mirroring the usual slowing down of thermally activated dipolar dynamics for low temperatures. Surprisingly, a non-canonical increase of peak frequency is observed when cooling further, which implies an acceleration of the involved molecular dynamics. Such behavior is typical for order–disorder ferroelectrics below the phase transition55.

It should be noted that the detected dielectric relaxation behavior (Fig. 2) in some respects resembles that of so-called relaxor ferroelectrics, which can be regarded as nano-ordered ferroelectrics with a smeared-out diffusive phase transition56. There are model simulations predicting typical signatures of relaxors—dispersive maxima due to polar nanoregions in the temperature-dependent dielectric permittivity—in bulk water and in proteins hydration shells57,58. Concerning the permittivity, the main distinction between a relaxor and an order–disorder ferroelectric phase transition observed in the present data is the critical behavior of the relaxation time τ (or critical slowing down of the peak frequency νp 1/τ; Fig. 3), which is typical for order–disorder ferroelectrics. In a relaxor, it should exhibit glass-like freezing following the Vogel–Fulcher–Tammann law instead59,60.

Fig. 3: Dielectric and thermodynamic indications of a phase transition at T0 = 3 K.
figure 3

Temperature dependences of peak frequency (dark blue circles) and dielectric strength ∆ε (purple circles) of the relaxational excitation (Fig. 2) observed in the hydrous cordierite crystal for the polarization E || a. The dark blue lines correspond to the fits according to ~ exp{−Ea/kBT}*(TT0) with Ea = 6.2 meV above T0 (T0 = 3 K; solid line) and ~ (T0T) below T0 (dashed line), as described in the text. The solid black triangles show the temperature dependence of the excess specific heat of water molecules hosted by cordierite. Δcp is calculated as difference between the specific heat of hydrated and dehydrated crystals. The open black triangles show the specific heat plotted as Δcp/T to stress the low-temperature part. Inset: temperature dependences of the peak frequency νp of the excitation and of the dipole relaxation time τ = (2πνp)−1 in a linear temperature scale. The error bars correspond to the ranges of the data that provide satisfactory description of the temperature-dependent relaxational excitation seen in the spectra of the complex permittivity.

Discussion

The temperature evolution of the detected relaxation process is quantitatively summarized in Fig. 3 by presenting the peak frequency νp and the dielectric strength ∆ε of absorption band. At 3 K < T < 8 K, the values of ∆ε were determined via fits of the RF spectra of Fig. 2b and c with the Havriliak–Negami function,

$$\varepsilon ^ \ast \left( \nu \right) = \varepsilon _\infty + {\Delta}\varepsilon \left[ {1 + \left( {{\rm{i}}\omega /\omega _{\mathrm{R}}} \right)^{1 - \alpha }} \right]^{ - \beta },$$
(2)

where ω is the angular frequency, ωR is the angular relaxation frequency and the exponents α and β characterize the broadness and asymmetry of the relaxational band, respectively. Fitting results for three selected temperatures are shown in Fig. 2b, c by short-dashed lines. In the lowest temperature range, 0.3 K < T0 < 3 K, the values of ∆ε were determined by phenomenologically describing the RF spectra of Fig.2a, b with the sum of two Havriliak–Negami relaxational terms with the dielectric contributions ∆ε1 and ∆ε2. Around T0 ≈ 3 K the peak frequency νp and the relaxation time τ = (2πνp)−1 (inset in Fig. 3) exhibit a pronounced minimum and maximum, respectively, which are clear signatures of an order–disorder ferroelectric phase transition55. The temperature dependence of the dielectric strength Δε exhibits a characteristic peak at the transition temperature, again typical for order–disorder ferroelectrics. In addition, the specific heat clearly reveals an anomaly around T0.

Since the electric dipole–dipole interaction energy of water molecules along the channels is Ud–d ≈ 22 meV47, at room temperature (kBT = 26 meV) the H2O molecules are mostly decoupled. The situation changes with cooling: the H2O dimers, trimers, etc., start to form in the channels47, leading to the Curie–Weiss behavior of Eq. (1) with C = 390 K and TC = −15 K (solid blue line in Fig.2a). The negative Curie temperature indicates predominantly antiferroelectric correlations between the dipoles in channels55,61, as indicated by our DFT-MD and Monte Carlo simulations, see below. Despite strong interactions, we do not observe any detectable signs of a phase transition associated with the appearance of long-range order of the H2O dipoles at elevated temperatures, most probably due to the incomplete filling (see Sample preparation and characterization section), as vacancies disrupt the one-dimensional water chains.

Below T ≈ 40 K, however, the thermal energy (kBT = 3.5 meV) becomes comparable to the in-plane dipolar coupling strength Ud–d ≈ 3 meV47; deviations from the Curie–Weiss dependence begin to evidence dipolar interactions within the ab-planes. At these temperatures the main contribution to the relaxational peak stems from complexes containing water dipoles coupled both along the channels c-axis and within the ab-planes. Due to polar ordering, the relaxational dynamics of these complexes freezes out at T0 ≈ 3 K leading to a strong decrease of the relaxational strength of the band in the ordered state. Concurrently the peak frequency of the band exhibits a V-shaped temperature dependence (Fig. 3). The softening of the peak frequency at 3 K < T < 10 K can be fitted by νp~ exp{−Ea/kBT}*(TT0) with Ea = 6.2 meV and T0 ≈ 3 K (blue solid line in Fig. 3). The Arrhenius factor corresponds to the non-interacting thermally activated relaxations and the (TT0) term describes the critical slowing down of the relaxational response of molecular complexes. At T < T0, the peak frequency follows the T0T dependence (blue dashed line in Fig. 3). In all details, the observed features correspond to the behavior of order–disorder ferroelectrics62,63,64,65.

The appearance of an ordered state of H2O molecular dipoles below T0 = 3 K is conclusively demonstrated by the pyrocurrent that peaks at the same temperature, and the polarization appearing at T < T0, see Fig. 4. The value of low-temperature polarization P ≈ 3 nC cm−2 allows us to estimate the number of dipoles that contribute to this polarization N = P/pa ≈ 1.4×1019 cm−3. (Here pa is component of H2O dipole along the a-axis). The obtained value is about two orders of magnitude smaller than the total concentration of water molecules in the crystal ≈1.9×1021 cm−3. The reason for this discrepancy lies in the fact that the H2O dipoles form ferroelectric domains at temperatures below T0 = 3 K but also acquire antiferroelectric order along the c-axis. An estimate of the field binding two collinear dipoles separated by a distance ≈ 10 Å within domain provides a value of 750 kV cm−147 that is comparable to crystalline internal fields. Then, the fields of 8 kV cm−1 used in our experiments will be able to polarize only small fraction of the dipoles (in the present case about 1%) that are located at the domains boundaries, close to the “defects” (empty cages) or isolated dipoles. We believe that for the same reasons we were not able to detect any meaningful signs of dielectric nonlinearity (hysteresis) in the fields up to 15 kV cm−1 because they are also not sufficient to affect the electric coupling between H2O dipoles.

Fig. 4: Pyrocurrent and polarization indicate a phase transition at T0 = 3 K.
figure 4

Temperature dependences of pyrocurrent (red) measured for E || a polarization while heating in zero electric field after cooling hydrous cordierite crystal in external electric field 8 kV cm−1. The temperature-dependent polarization (blue) was calculated from pyrocurrent. Open dots correspond to the pyrocurrent measured during heating in zero electric field after cooling the crystal in zero field.

To provide a more detailed picture of the dipolar ordering, we performed density functional theory molecular dynamics simulations of two H2O molecules located next to each other along the channel c-axis and within the ab-plane. The simulations were performed for various temperatures and yield averaged snapshots of the oxygen and two protons of the corresponding two molecules in 1 fs steps during a period of 15 ps (see Supplementary Movie 16). Figure 5 displays the ab-projection of the positions of oxygen and two protons in two adjacent cages along the channel and in neighboring channels, respectively. At room temperature, hardly any correlations between the sites can be identified, the molecules are nearly independent. Upon cooling, however, the disorder due to thermal motion weakens and the molecules become progressively more confined. The important observation is that the molecules prefer to orient their dipole moments antiparallel within the channels while the dipoles align parallel within the ab-planes. Analogous simulations performed for the larger system of four adjacent H2O molecules in a channel, in plane and occupying all four cages of cordierite crystal unit cell, see Supplementary Figs. 24, do not qualitatively change the preferred configurations of dipoles. Only the potentials experienced by the molecules acquire a slightly more complex shape with local minima for oxygen and protons. The results of our Monte Carlo analysis of the dynamics of N = 3072 interacting water molecules (see Monte Carlo simulations methods section) are in full agreement with the conclusions drawn from the DFT-MD simulations. As demonstrated in Fig. 6, the additional finding from the Monte Carlo analysis is that in order to minimize the energy of dipolar system, the molecules tend to form in-plane ferroelectric domains with a polarization predominantly along the b-axis and alternating its sign along the c-direction. We calculated the distribution of domains sizes as a function of their lengths at temperature T = 0.001 K (see Supplementary Fig. 5). Domains are considered as arrays of collinear (along the b-axis) dipoles with their boundaries given by dipoles of different directions, by defects (empty cages) or by boundaries of the sample. The obtained mean size of the domains is 1.75 and 2.28 lattice distances (see Supplementary Fig. 5) along the a-axis and b-axis, respectively; after multiplication by the translation lattice vectors the average domains size is 1.496 nm×2.220 nm = 3.32 nm2. Figure 7 displays the temperature dependences of the antiferroelectric order parameter along the a-axis and a-axis dielectric susceptibility, both clearly demonstrating a phase transition at T0 = 3 K. In other words, our DFT-MD and Monte Carlo results confirm a three-dimensional low-temperature order of the H2O molecular network. The water molecules in cordierite form a dipolar lattice with ferroelectric domains in ab-planes staggered antiferroelectrically along the channels c-axis, as shown in Fig. 8.

Fig. 5: Temperature-dependent dynamics of nanoconfined water molecules.
figure 5

Density functional theory molecular dynamics simulation during 15 ps and at five temperatures of two water molecules within two nanocages of cordierite crystal lattice, located next to each other along the channel c-axis (left) and within the ab-plane (right). The panels display projections of oxygen ions (O) and two protons (H1, H2) of the H2O molecules on the ab-plane; the c-axis is perpendicular to the figure plane. Red arrows represent the dipole moments of the molecules that show antiferroelectric alignment for molecules located along the channel and ferroelectric alignment for molecules in the ab-plane. Numbers correspond to the coordinates of atoms within the cages in fractions of translation vectors of cordierite unit cell.

Fig. 6: Configuration of water molecular dipoles at different temperatures.
figure 6

We plot the configurations of water molecular dipoles in ab and ac cross sections at different temperatures obtained from Monte Carlo simulations. Red and blue colors correspond to directions parallel (blue) and antiparallel (red) to the b-axis; bright and faded tones correspond to directions along (bright) and against (faded) the a-axis. Since the dipoles lie within the ab-planes and have zero component along the c-axis, arrows in the ac-panels are artificially turned to indicate their directions in the ab-planes. For the analysis, angles of ±20° relative to the b-axis were taken from DFT-MD simulations. Filling factor is 75%.

Fig. 7: Dielectric susceptibility and antiferroelectric order parameter.
figure 7

Temperature dependences of the dielectric susceptibility (red) and antiferroelectric order parameter (blue) of nanoconfined H2O molecules along the a-axis. The dielectric susceptibility is defined through the polarization \(P_{\mathrm{a}} = \left( {1/V} \right)\mathop {\sum}\nolimits_{i,j,k} {p_{ijk}^{\mathrm{a}}}\), where V is the crystal volume, the sum is taken over all lattice points (i, j, k) corresponding to the lattice axes (a, b, c), \(p_{ijk}^{\mathrm{a}}\) is the component of the water dipole p along the a-axis at a position (i, j, k), the electric field component Ea along the a-axis enters the corresponding expression for dielectric susceptibility via \(\chi _{{\mathrm{aa}}} = P_{\mathrm{a}}/\varepsilon _0E_{\mathrm{a}}\), ε0 is the permittivity of vacuum. The antiferroelectric order parameter is calculated as \(\eta _{AF}^{\mathrm{a}} = \left( {1/N\left| {p^{\mathrm{a}}} \right|} \right)\mathop {\sum}\nolimits_{i,j,k} {\left( { - 1} \right)^kp_{ijk}^{\mathrm{a}}}\), where N is the number of dipoles in the crystal. The dependences were obtained by Monte Carlo simulations for the electric-dipolar system of water molecules in hydrous cordierite crystal with N = 3072 and filling factor 75%.

Fig. 8: Configuration of H2O dipole moments from Monte Carlo simulations.
figure 8

We show the configuration of the molecular dipoles of water in cordierite crystalline matrix at zero temperature with two filling factors of water molecules. Red and blue colors correspond to directions along (blue) and against (red) the b-axis. a 100% filling factor. Dipoles are oriented ferroelectrically in the ab-planes and antiferroelectrically along the c-axis. b 75% filling factor corresponding to the studied hydrous cordierite crystal.

Based on our comprehensive experimental and numerical investigations of hydrous cordierite crystals we firmly conclude that H2O molecules form a highly ordered dipolar lattice at low temperatures. When confining water molecules to the anisotropic nanosized channels composed by the ionic lattice of the crystal, they undergo a ferroelectric order–disorder type phase transition near T0 = 3 K with characteristic features in the temperature dependence of dielectric permittivity, specific heat, pyroelectric current and polarization. Density functional theory molecular dynamics and Monte Carlo simulations indicate the antiferroelectric coupling of H2O molecular dipoles along the nanochannel c-axis and their ferroelectric coupling within the ab-planes. As a result, the low-temperature phase is characterized by in-plane ferroelectric domains ordered antiferroelectrically along the nanochannels direction. This long-sought polar phase transition in a system of coupled dipolar water molecules demonstrates that hydrous crystals provide an ideal workbench for studies of different phases and phase transitions in the electric counterpart of magnetic spin systems.

Methods

Sample preparation and characterization

Natural cordierite crystal from India (the detailed location is unknown) have been carefully selected and studied. According to visual inspection with an optical microscope at low magnification (×10), the crystal was free of impurities or foreign phases. 10 independent measurements via microprobe analysis (JEOL JXA-8100, Analytical Center for multi-elemental and isotope research SB RAS) provided with the following chemical composition in mass %: Na2O—0.204; MgO—13.350; MnO—0.017; K2O—0.008; CaO—0.021; SiO2—49.296; Al2O3—33.140; FeO—1.584; Loss on ignition—2.200 (total—99.821). Recalculation of this chemical analysis on the basis of 18 oxygen atoms per formula gives the following chemical formula:

(Mg1.643Mn0.001Fe0.137Al0.104Ca0.002)Σ = 1.887Al3 (Al0.921Si5.079)Σ = 6O18(H2O)0.756 (Na0.041K0.001)Σ = 0.042.

This recalculation assumes that the total weight loss on ignition is due to water loss, which is generally not the case, since some unknown amount of CO2 is also presents in cordierite. Thus, the water content specified in the formula is an estimate from above. From the formula, the filling factors of water-I and water-II47 can be estimated as ~75.6% and ~4.2%, respectively.

To solve the complex problem of identification of water molecules in the studied cordierite crystal, an accurate X-ray diffraction analysis was carried out. The sets of diffraction reflection intensities at different temperatures were collected on a Rigaku Oxford Diffraction Xcalibur diffractometer equipped with an EOS S2 detector. Our refinement of the cordierite structure provided with the following results. Space group is Cccm, Z = 4. Maximum angle of scattering of reflections is θ = 74.3°. Unit cell parameters are a = 17.1004(3) Å, b = 9.7399(2) Å, c = 9.3268 (2) Å. The residuals that describe the differences between the experimental and modeling data are R/wR = 1.16/1.35%. Extremes of the difference Fourier synthesis of electron density are Δρminρmax = –0.17/+0.29 electrons per Å3 for 5693 symmetry-independent reflections. Structural results are presented in Supplementary Table 1. It was found that water molecules are located in the region (0, 0, ¼). The vector connecting protons of the molecule is almost parallel to the crystallographic c-axis, and the vector of the molecular dipole moment is slightly tilted relative to the b-axis, in agreement with the present DFT-MD results.

Experimental techniques

From the crystal, several slices were cut with the a, b, and c crystallographic axes lying within or perpendicular to their planes; these geometries allowed to measure the broad-band dielectric response of the samples in all three principal polarizations of the electric field vector E of the probing electromagnetic radiation. For the measurements at frequencies from 1 Hz up to the microwave range, we used a Novocontrol Alpha AN High Performance Frequency Analyzer equipped with a He-flow cryostat JANIS ST-100, an Andeen-Hagerling 2500A capacitance bridge connected to the cryostat utilizing a single-shot 3He insert with embedded coaxial cables, and a coaxial reflectometric technique employing an impedance analyzer Keysight 4991B66. A time-domain TeraView 3000 terahertz spectrometer was employed to determine the dielectric spectra up to few terahertz. Dielectric hysteresis loops were studied at frequencies 1–50 Hz using the standard Sawyer-Tower Bridge. The pyrocurrent measurements (0.5–2 K min−1 heating after cooling with field up to 8 kV cm−1) were realized using the Electrometer High Resistance Meter KEITHLEY 6517B. The polarization was obtained from the pyrocurrent data. The dielectric experiments were complemented by measurements of the heat capacity in the relaxation method employing a PPMS system (Quantum Design). In all experiments, measurements on dehydrated samples (annealing at a temperature 1000 °C for 8 h according to ref. 67) allowed us to extract the characteristics determined exclusively by water molecules. Complete dehydration was monitored by the weight loss of the specimen and by disappearance of internal vibrations of the H2O molecule in the infrared spectra.

Density functional theory molecular dynamics simulation

Ab initio DFT calculations were carried out using the Vienna Ab initio Simulation Package (VASP)68,69 with plane wave basis sets cut off at 400 eV, PAW pseudopotentials70, and Perdew–Burke–Ernzerhof (PBE) exchange-correlation functional71. The Brillouin zone integration was carried out using a 1x2x2 k-point sampling and Gaussian smearing. One unit cell of cordierite of 116 atoms was simulated with initial geometry taken from X-ray diffraction (Supplementary Table 1). For all calculations lattice parameters were fixed to the experimentally observed values (see section Sample preparation and characterization above) and periodic boundary conditions were applied. One unit cell of cordierite forms four nanocages, and several types of nanocage filling were simulated. At the first stage, one water molecule and two water molecules within nanocages arranged along the c-axis and within ab-plane were studied. For each geometry, molecular dynamics simulations were performed at 10 K, 50 K, 75 K, 100 K, 150 K, 300 K for 15 ps or 20 ps with 1 fs time step using Nose–Hoover thermostat. Additionally, molecular analysis was performed for four molecules occupying all four cages of the unit cell and for large system of 1x1x2 and 1x2x1 super cell with four water molecules arranged next to each other along the channel c-axis and within the ab-plane, respectively. In each procedure, the first 1.5 ps of simulation were considered as thermalization and removed from the subsequent analysis. To find out whether van der Waals (vdW) forces play any significant role in interaction between the water molecules and the ions of nanocage, we performed test calculations using the non-local dispersion corrected vdWDF2 functional72 (Supplementary Fig. 6, Supplementary Note 1).

Monte Carlo simulations

Metropolis algorithm was used to simulate at different temperatures the behavior of the dipole system governed by the dipole–dipole interaction Hamiltonian \(H = \left( {8{\uppi}\varepsilon _{\mathrm{r}}\varepsilon _0} \right)^{ - 1}\mathop {\sum}\nolimits_{ij} {r_{ij}^{ - 3}\left( {{\mathbf{p}}_i{\mathbf{p}}_j - 3\left( {{\mathbf{p}}_i{\mathbf{n}}_{ij}} \right)\left( {{\mathbf{p}}_j{\mathbf{n}}_{ij}} \right)} \right)}\). Here, ε0 is the dielectric constant, εr = 5 is the relative dielectric permittivity due to other degrees of freedom than the water dipoles, rij is the distance between two dipoles pi and pj, and nij is the unit vector between them. The presented results were obtained for a sample with 16 lattice sites along each axis (a total of 4096 sites) and a 75% dipole filling factor with N = 3072 dipoles in total (25% of sites were free of dipoles, we called them defects). Free boundary conditions were applied. These results were well reproduced in simulations with a finite interaction radius cut-off rc (with up to the 4th next-to-nearest neighbors in the ab-plane, which corresponds to 640 lattice sites in the interaction sphere) with periodic boundary conditions. However, in the latter case, the dipoles configuration at the lowest temperature usually consists of domains of radius rc even in the absence of defects and converges slowly to the ground state. Therefore, we avoided this method in order to see the dipoles configurations in the presence of defects more reliably. The dielectric susceptibility along each axis α was calculated from fluctuations of the average dipole \(p_\alpha = N^{ - 1}\mathop {\sum}\nolimits_i^N {p_{\alpha i}}\) as \(\chi _\alpha = Np_0^2\left( {\upsilon _0\varepsilon _0k_BT} \right)^{ - 1}\left( {\left\langle {p_\alpha ^2} \right\rangle - \left\langle {p_\alpha } \right\rangle ^2} \right)\). Here, N is the number of dipoles in a simulation, p0 = 1.85 D is the water molecular dipole, υ0 is the volume per dipole, kB is the Boltzmann constant. For each temperature, the number of Monte Carlo steps per spin was 3500 (the first 500 were attributed to thermalization). The results were averaged over 30 samples with different randomly generated defect configurations. From DFT-MD analysis it was deduced that at low temperatures the dipoles tend to orient at certain angle relative to the b-axis. An angle of 20° was taken for all simulations.