Dielectric ordering of water molecules arranged in a dipolar lattice

Intermolecular hydrogen bonds impede long-range (anti-)ferroelectric order of water. We confine H2O molecules in nanosized cages formed by ions of a dielectric crystal. Arranging them in channels at a distance of ~5 Å with an interchannel separation of ~10 Å prevents the formation of hydrogen networks while electric dipole-dipole interactions remain effective. Here, we present measurements of the temperature-dependent dielectric permittivity, pyrocurrent, electric polarization and specific heat that indicate an order-disorder ferroelectric phase transition at T0 ≈ 3 K in the water dipolar lattice. Ab initio molecular dynamics and classical Monte Carlo simulations reveal that at low temperatures the water molecules form ferroelectric domains in the ab-plane that order antiferroelectrically along the channel direction. This way we achieve the long-standing goal of arranging water molecules in polar order. This is not only of high relevance in various natural systems but might open an avenue towards future applications in biocompatible nanoelectronics.

F or decades intense research activities tackle the question whether water molecules with their rather strong dipole moment of p 0 = 1.85 Debye can condense into a ferroelectrically or antiferroelectrically ordered state. In pure liquid water or in H 2 O 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 rules [1][2][3] . Proton ordered phases of water ice can be induced by defects introduced via doping with impurities or via irradiation 3,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 H 2 O 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 H 2 O molecules are found in fully or partially confined states in cells, membrane channels, and proteins hydration shells [5][6][7][8] , which can influence the functioning of living organisms substantially.
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 silicates 42,43 , the special appeal of these crystals is that they contain just separate H 2 O 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 moments 44,45 . In addition, introducing dopants like Na and K ions at bottlenecks that separate nanovoids can lead to a polarization of H 2 O molecules and the formation of statically polarized so-called water-II molecules. Having also the possibility to encapsulate heavy water molecules (D 2 O, DHO), one gets a unique "laboratory" for studies of various excitations, phases and phase transitions in the water dipolar network. Very recently it was reported 46-51 that H 2 O 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 lattice 52 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 H 2 O 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 T 0 ≈ 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) 2 Al 4 Si 5 O 18 is ideally suited for these requirements as it is structurally similar to beryl, except for the lack of rotational symmetry. Beryl contains Si 6 O 18 rings, while in cordierite two of six silicon atoms are replaced by aluminum, leading to (Si,Al) 6 O 18 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 Cccm 53 . There are indications that the H 2 O dipole moments in cordierite are preferably oriented parallel to the baxis 54 . 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.
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 T 0 ≈ 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 H 2 O 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 H 2 O 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: Here, ε ∞ is the high-frequency dielectric constant, C is the Curie constant and T C is the Curie temperature. The lowtemperature 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 T 0 ≈ 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 transition 55 .
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 transition 56 . 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 shells 57,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 instead 59,60 .

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, 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 < T 0 < 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 T 0 ≈ 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 transition 55 . 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 T 0 .
Since the electric dipole-dipole interaction energy of water molecules along the channels is U d-d ≈ 22 meV 47 , at room temperature (k B T = 26 meV) the H 2 O molecules are mostly decoupled. The situation changes with cooling: the H 2 O dimers, trimers, etc., start to form in the channels 47  Below T ≈ 40 K, however, the thermal energy (k B T = 3.5 meV) becomes comparable to the in-plane dipolar coupling strength U d-d ≈ 3 meV 47 ; 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 T 0 ≈ 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~e xp{−E a /k B T}*(T−T 0 ) with E a = 6.2 meV and T 0 ≈ 3 K (blue solid line in Fig. 3). The Arrhenius factor corresponds to the non-interacting thermally activated relaxations and the (T-T 0 ) term describes the critical slowing down of the relaxational response of molecular complexes. At T < T 0 , the peak frequency follows the T 0 -T dependence (blue dashed line in Fig. 3). In all details, the observed features correspond to the behavior of order-disorder ferroelectrics [62][63][64][65] .
The appearance of an ordered state of H 2 O molecular dipoles below T 0 = 3 K is conclusively demonstrated by the pyrocurrent that peaks at the same temperature, and the polarization appearing at T < T 0 , 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/p a ≈ 1.4×10 19 cm −3 . (Here p a is component of H 2 O 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×10 21 cm −3 . The reason for this discrepancy lies in the fact that the H 2 O dipoles form ferroelectric domains at temperatures below T 0 = 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 −1 47 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 H 2 O dipoles.
To provide a more detailed picture of the dipolar ordering, we performed density functional theory molecular dynamics simulations of two H 2 O 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   (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. The dark blue lines correspond to the fits according to~exp{−E a /k B T}*(T −T 0 ) with E a = 6.2 meV above T 0 (T 0 = 3 K; solid line) and~(T 0 −T) below T 0 (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. Δc p is calculated as difference between the specific heat of hydrated and dehydrated crystals. The open black triangles show the specific heat plotted as Δc p /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.
of the oxygen and two protons of the corresponding two molecules in 1 fs steps during a period of 15 ps (see Supplementary Movie 1-6). 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 H 2 O molecules in a channel, in plane and occupying all four cages of cordierite crystal unit cell, see Supplementary Figs. 2-4, 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 nm 2 . 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 T 0 = 3 K. In other words, our DFT-MD and Monte Carlo results confirm a three-dimensional low-temperature order of the H 2 O 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.
Based on our comprehensive experimental and numerical investigations of hydrous cordierite crystals we firmly conclude that H 2 O 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 T 0 = 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 H 2 O 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. 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 CO 2 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-II 47 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 singleshot 3He insert with embedded coaxial cables, and a coaxial reflectometric technique employing an impedance analyzer Keysight 4991B 66 . 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  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 pseudopotentials 70 , and Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional 71 . 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 functional 72 ( 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 Here, ε 0 is the dielectric constant, ε r = 5 is the relative dielectric permittivity due to other degrees of freedom than the water dipoles, r ij is the distance between two dipoles p i and p j , and n ij 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 r c (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 r c 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 α ¼ N À1 P N i p αi as Here, N is the number of dipoles in a simulation, p 0 = 1.85 D is the water molecular dipole, υ 0 is the volume per dipole, k B 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.  Temperature dependences of the dielectric susceptibility (red) and antiferroelectric order parameter (blue) of nanoconfined H 2 O molecules along the a-axis. The dielectric susceptibility is defined through the polarization P a ¼ 1=V À ÁP i;j;k p a ijk , 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 a ijk is the component of the water dipole p along the a-axis at a position (i, j, k), the electric field component E a along the a-axis enters the corresponding expression for dielectric susceptibility via χ aa ¼ P a =ε 0 E a , ε 0 is the permittivity of vacuum. The antiferroelectric order parameter is calculated as η a AF ¼ 1=N p a j j À ÁP i;j;k À1 ð Þ k p a ijk , 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%. 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.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.