Proton Ordering of Cubic Ice Ic: Spectroscopy and Computer Simulations.

Several proton-disordered crystalline ice structures are known to proton order at sufficiently low temperatures, provided that the right preparation procedure is used. For cubic ice, ice Ic, however, no proton ordering has been observed so far. Here, we subject ice Ic to an experimental protocol similar to that used to proton order hexagonal ice. In situ FT-IR spectroscopy carried out during this procedure reveals that the librational band of the spectrum narrows and acquires a structure that is observed neither in proton-disordered ice Ic nor in ice XI, the proton-ordered variant of hexagonal ice. On the basis of vibrational spectra computed for ice Ic and four of its proton-ordered variants using classical molecular dynamics and ab initio simulations, we conclude that the features of our experimental spectra are due to partial proton ordering, providing the first evidence of proton ordering in cubic ice. We further find that the proton-ordered structure with the lowest energy is ferroelectric, while the structure with the second lowest energy is weakly ferroelectric. Both structures fit the experimental spectral similarly well such that no unique assignment of proton order is possible based on our results.

Introduction

Despite its simple molecular structure, water has a remarkably complex phase diagram. Application of pressure produces a variety of different solid ice phases with densities considerably higher than that of ordinary hexagonal ice, ice Ih. To date, 16 thermodynamically stable or metastable crystalline phases (labeled with Roman numerals as Ih, Ic, II, III, ..., XV)[1−3] and several amorphous phases[4−6] have been discovered. While in some crystalline ice structures, including ice Ih as well as cubic ice, ice Ic, only the oxygen atoms form a regular lattice and the protons are disordered, in other ice phases, such as ice II, also the protons are arranged in a regular way. Indeed, there exist pairs of ice structures, such as ice Ih and its protonically ordered counterpart ice XI, which have nearly identical oxygen sublattices but differ in their proton order. Recently, Salzmann and collaborators have identified experimentally[2,3] several previously unknown ice phases, completing the pairings V/XIII, XII/XIV, and VI/XV, which differ only in proton ordering.

Proton ordering transitions in crystalline ices take place at rather low temperatures and are severely hampered by slow transformation kinetics. The highest ordering transition temperatures have been found for ice VII and ice III, which order below T ≈ 270 K at p > 2 GPa to ice VIII[7] and below T ≈ 170 K at p ≈ 0.3 GPa to ice IX,[8,9] respectively. These two hydrogen ordering transitions in ice are the only ones known to take place in the absence of a catalyst. Other ordering transitions in the intermediate pressure regime at p ≈ 0.5–1.5 GPa take place at about T = 100–150 K, namely, ice V ⇌ ice XIII,[2] ice VI ⇌ ice XV,[3] and ice XII ⇌ ice XIV.[2] These three transitions take place at an observable rate by using HCl as a dopant, which increases the reorientation rates of individual water molecules in the ice lattice, presumably by producing rotational Bjerrum L-defects.[1] Ambient-pressure hexagonal ice, ice Ih, shows the lowest experimentally found ordering transition temperature of T ≈ 72 K.[10−17] Such proton-ordered ice XI is typically produced from ice Ih using hydroxide doping, for example, by freezing a 0.1 M KOH solution.

Because ice Ih and Ic have identical first coordination shells and differ only in the stacking of the planes orthogonal to the c-axis, one would expect that ice Ih and Ic proton order under similar conditions. However, the proton-ordered counterpart of cubic ice Ic has not been observed to date. To study the energetics of proton arrangements in ice Ic, Lekner computed the electrostatic energies of all 90 proton configurations satisfying the ice rules in a periodically replicated unit cell of 8 water molecules.[18,19] Due to the degeneracy of the Coulomb energy, there are only four classes of configurations with different energies, examples of which are shown in Figure 1. Out of these, the perfectly ordered antiferroelectric proton configuration has the lowest electrostatic energy, while the ferroelectric structure has the highest electrostatic energy. Weakly ferroelectric configurations have intermediate energies.

Figure 1

Four distinct proton arrangements in a unit cell of cubic ice containing eight water molecules[18] with hydrogen bonds (in green). Each configuration is a representative of one of the four symmetry-inequivalent proton arrangements with different Coulomb energies that exist in a unit cell of this size. Configuration a (space group I41md) is ferroelectric, configurations b (Pna21) and c (P41) are weakly ferroelectric, and configuration d (P41212) is antiferroelectric.

This general trend is confirmed by calculations for larger unit cells such that, on the basis of these results, thermodynamically one expects a transition to an antiferroelectric phase for ice Ic at sufficiently low temperatures. The purely electrostatic calculations of Lekner, however, are contrasted by recent ab initio (density functional theory) calculations,[20] according to which the ferroelectrically ordered structure is the energetically most stable one. Similar differences between ab initio results and calculations based on empirical force fields have been reported also for the energetic ordering of various proton-ordered forms of hexagonal ice.[21] Most likely, such discrepancies are due to polarization effects,[22] which are fully taken into account in ab initio calculations but are almost entirely neglected by most empirical potentials.

In this paper, we investigate the proton ordering of ice Ic using a combination of spectroscopic experiments and computer simulations. As discussed in detail below, the shape of the librational band of the IR spectra measured in our work, interpreted on the basis of computed IR spectra, indicates that hydroxide-doped ice Ic partially proton orders following an experimental protocol similar to that used previously to proton order ice Ih. Interestingly, the combination of experiment and theory suggests that cubic ice indeed shows a tendency to ferroelectric proton order; however, kinetics may strongly influence the ordering transition.

The remainder of this article is organized as follows. In section 2, we explain the computational model and methods used to carry out the simulations and report their results. Experimental procedures and results are described in section 3 and discussed and interpreted in section 4 based on the spectra computed in our simulations.

Simulations

In order to detect proton-ordered cubic ice in the laboratory, it is important to be capable of distinguishing the various degrees of protonic order. One possible way to do that is via the signature of proton order on the vibrational spectrum obtained in infrared (IR) or Raman spectroscopy experiments, methods that have been employed before to study proton order in ice.[23−25] In the present work, we used molecular dynamics (MD) simulations complemented with ab initio calculations to study the structural and dynamical properties of various proton-ordered ice Ic polymorphs. In particular, we determined radial distribution functions and IR spectra for different proton-ordered cubic ice candidates and compared them to the ones for hexagonal ice Ih and its proton-ordered counterpart, ice XI. In doing so, we focused on translational and librational modes as these modes carry the most information about proton ordering and provide the basis for the interpretation of the experimental IR spectra presented in section 3.

Model and Methods

All of our MD simulations are performed using the TIP4P/ice model,[26] which reproduces all ice phases consisting of intact water molecules and leads to a phase diagram with the correct topology and coexistence lines that are only slightly displaced with respect to the experimental phase diagram.[27] Because in the TIP4P/ice model water molecules are rigid, only motions that do not involve intramolecular vibrations can be studied such that our IR spectra are limited to the low-frequency bands corresponding to translational and librational modes.

We integrated the equations of motion with a quaternion-based integrator that maintains the rigid geometry of the water molecules. In particular, we carried out MD simulations in the NVT and isotropic NpT ensembles using a slightly modified version of the Verlet-like integrator proposed by Kamberaj et al.,[28] based on the Trotter decomposition schemes applied by Miller et al.[29] and Martyna et al.[30] In this algorithm, the canonical and isothermal–isobaric MD is implemented through thermostat chains based on the Nosé–Hoover[31,32] and the Andersen[33] approaches. Electrostatic interactions were treated with Ewald summation.[34]

The initial configurations of the different proton-ordered phases were generated by periodically replicating the unit cells given in ref (18) for cubic ice and in refs (21) and (12) for ice XI, respectively. The proton-disordered counterparts were set up following the procedure suggested by Rahman and Stillinger,[35] which starts from a hexagonal or cubic crystal with perfect proton order and protons located on the connecting lines between neighboring oxygens. The proton order is then disrupted by shifting protons from a position close to one oxygen to the position close to the neighboring oxygen originally accepting the hydrogen bond involving the shifted hydrogen. In order to keep the water molecules in the sample intact, this shifting operation has to be carried out along closed loops of hydrogen bonds. This sequence of basic steps is repeated until every water molecule has been touched several times and a vanishing total dipole moment is achieved. At the end of the procedure, the ideal water molecules, which have a geometry consistent with perfect tetrahedral coordination, are replaced by molecules with TIP4P geometry.

To quantify the structural and dynamical diversity of various degrees of proton order, we have calculated energies, pair correlation functions, as well as IR absorption spectra. In the present case, we consider the pair correlation functions gOO(r), gHH(r), and gOH(r) between pairs of oxygen atoms, pairs of hydrogen atoms, and oxygen and hydrogen atoms, respectively. All of these functions can be extracted from data obtained in neutron diffraction experiments.[36] IR absorption spectra are calculated in the classical approximation[37,38] as the Fourier transform of the time autocorrelation function ⟨M(0)·M(t)⟩ of the total dipole moment Mwhere ω is the vibrational frequency and angular brackets ⟨···⟩ indicate a time or ensemble average. The dipole–dipole correlation function can be written in terms of the dipole moments μ of individual water molecules, which are determined by the magnitude, sign, and location of the charges on the TIP4P/ice molecules as well as by their orientation, ⟨M(0)·M(t)⟩ = ∑ ⟨μ(0)·μ(t)⟩. Note that while the cross correlation terms ⟨μ(0)·μ(t)⟩, i ≠ j, are negligible for high-frequency modes associated with intramolecular stretching and bending motions, they play an important role for the form of the spectrum in the range characteristic for translations and librations involving collective motions of multiple molecules.[39,40]

We complemented our simulations based on the TIP4P/ice model with energies and spectra computed ab initio using the Vienna ab initio simulation package (VASP) and PAW[41] potentials in the implementation of Kresse and Joubert.[42] The outermost core radii for the O and H potentials are 1.52 and 1.1 au, respectively (corresponding to the standard potentials distributed with the VASP package). All calculations were performed using the Perdew–Burke–Ernzerhof (PBE) functional,[43] as well as using van der Waals density functional theory (vdW-DFT).[44] We specifically used the vdW-DFT of Klimeš et al. termed “optPBE” (optimized PBE).[45,46] The Brillouin zone was sampled at 6 × 6 × 6 k-points. To determine the equilibrium volume of each structure, all internal parameters (including relative lattice parameters) were optimized at seven volumes around the equilibrium volume, and the energy versus volume curve was fitted using an equation of state. The vibrational frequencies were evaluated using finite differences; all symmetry-inequivalent atoms were displaced along symmetry-inequivalent directions, and the interatomic force constant matrix was completed using symmetry considerations.[47] The vibrational frequencies were determined at the Γ-point by diagonalization of the force constant matrix. The dipole activity was calculated from Born effective charge tensors for oxygen and hydrogen. For consistency with the MD simulations and to account for thermal expansion, the vibrational frequencies were evaluated at the average densities computed in the MD simulations. The vdW–DFT equilibrium volumes (at T = 0 K) are, however, only 3–4% smaller than the average MD volumes.

Simulation Results

To investigate the structure of ice Ih and ice Ic, we calculated radial distribution functions using NpT MD simulations of N = 1000 (for the parental cubic ice) and 896 (for the parental hexagonal ice) water molecules with orthogonal simulation boxes at temperature T = 170 K and pressure p = 0.1 bar. The proton-ordered counterparts were calculated using the same simulation boxes but are of lower space group symmetry, for example, orthorhombic symmetry in the case of ice XI. Similarly, proton ordering of cubic ice Ic leads to configurations that are not cubic but display a lower symmetry. When we speak about cubic and hexagonal ordered structures in the following, we always refer to the symmetry of the disordered parent structure. For the proton-disordered configurations, we also carried out simulations at T = 70 K, the temperature to which ice Ic is cooled in the experiments. Note that at a temperature of T = 170 K, initially proton-ordered configurations immediately disorder on the experimental time scale. In the simulations, however, no proton disordering is observed during the entire simulation because the simulation time is much shorter than the time scale of proton disordering. In all simulations, a time step of 2 fs and a Lennard-Jones cutoff of 3σ were used. The same cutoff was used for the real space part of the Ewald summation, where for the reciprocal space 1152 k vectors were employed. All simulations were performed for a hydrogen mass of mH = 1 and an oxygen mass of mO = 16. The total length of the simulations was 5 ns in each case.

Average densities and energies computed in these MD simulations are listed in Table 1 alongside the electrostatic energies computed by Lekner[18,19] as well as the energies determined in our ab initio simulations. The electrostatic calculations for the idealized structures and the MD simulations carried out at T = 170 and 70 K agree in the energetic ordering of the various structures, both finding that the antiferroelectric structure (d) has the lowest energy and the ferroelectric structure (a) has the highest energy while structures (b) and (c) have intermediate energies. The energies obtained using density functional theory, however, display the reverse energetic order but agree well with energies computed previously using DFT methods.[20] To double check the density functional theory data, we performed more accurate calculations using exact exchange and the random phase approximation (EXX-RPA) for the correlation energy.[48,49] These calculations were also extended to ice II, VIII, for which highly accurate diffusion Monte Carlo results are available.[50] In these cases (ice II, ice VIII, ice Ih), we found excellent agreement with the published diffusion Monte Carlo, validating the EXX-RPA.[51] The EXX-RPA results confirm the order of the cubic phases predicted by density functional theory; the ferroelectric structure (a) has the lowest energy, the antiferroelectric structure (d) has the highest energy, and the structures (b) and (c) have intermediate energies. This discrepancy, observed before for hexagonal ice,[21] is probably due to the neglect of polarization effects in the TIP4P/ice, which are expected to be particularly pronounced in ferroelectric structures, leading to a lowering of their energy. We would like to emphasize that predicting the correct energetic ordering is known to be difficult in the case of proton-ordered ices. Vega et al.[52] have shown that SPC/E, TIP4P, and TIP5P predict a transformation from ice Ih to the antiferroelectric Pna21 structure[53] below 70 K. The ferroelectric Cmc21 structure, which is obtained in experiments, is predicted to be the lowest-lying H-bond isomer only in the model by Nada and van der Eerden (NvdE).[54] It is also worth noting that the free energetic ordering of different proton-ordered configurations may also be significantly affected by higher-order multipoles, which are not taken into account in our TIP4P/ice simulations.[55]

Table 1

Average Densities ⟨ρ⟩ = mHN/⟨V⟩ and Average Potential Energies Per Molecule ⟨Epot⟩ for the Cubic and Hexagonal Ice Phases at Temperatures T = 70 and 170 K, respectively, and Pressure p = 0.1 bar along with the Coulomb Energies ECoul Calculated by Lekner[18,19] and the Cohesive Energies Per Molecule Computed Using Density Functional Theory (Ecoh) and Exact Exchange with the Random Phase Approximation (ERPA) at T = 0 and p = 0a

ice phase	⟨ρ⟩70K (kg/m3)	⟨ρ⟩170 K (kg/m3)	⟨Epot⟩70K (kJ/mol)	⟨Epot⟩170 K (kJ/mol)	ECoul (kJ/mol)	Ecoh (kJ/mol)	ERPA (kJ/mol)
ice Ih	930.17 (2)	919.00 (3)	0.0098 (2)	0.0169 (5)
ice XI (Cmc21)	930.31 (3)	919.12 (2)	0.0471 (2)	0.0547 (5)	0.750	–0.492
ice XI (Pna21)	929.95 (2)	919.25 (3)	–0.0195 (2)	–0.0250 (5)	–0.032	–0.096
ice Ic	932.11 (2)	919.89 (2)	0.0449 (2)	0.0524 (4)
ice Ic (ord. a, I41md)	929.32 (2)	918.39 (2)	0.2319 (2)	0.2299 (5)	0.816	–0.531	–0.521
ice Ic (ord. b, Pna21)	932.29 (2)	920.66 (3)	0.0364 (2)	0.0333 (4)	0.408	–0.261	–0.212
ice Ic (ord. c, P41)	932.24 (2)	920.51 (2)	0.0277 (1)	0.0261 (5)	0.204	–0.106	–0.068
ice Ic (ord. d, P41212)	932.57 (3)	920.63 (3)	0.0000	0.0000	0.000	0.000	0.000

All energies are stated with respect to the corresponding energy of the structure ice Ic (ord. d). For ⟨Epot⟩70K and ⟨Epot⟩170K, the total average potential energy per molecule of the reference structure (d) at T = 70 K is ⟨Epot⟩70K = −67.0872 (2) kJ/mol, and at T = 170 K ⟨Epot⟩170K = −64.3298 (3) kJ/mol. For Ecoh, the cohesive energy per molecule of the reference structure (d) with respect to infinitely separated water molecules is Ecoh = −63.813 kJ/mol. The statistical errors were calculated by block average analysis[34] and are given as single numbers in parentheses, which correspond to the error in the last digit.

The partial pair correlation functions of proton-disordered hexagonal ice, ice Ih, and its proton-ordered counterpart, ice XI, as well as disordered cubic ice, ice Ic, and its four proton-ordered candidates are given in Figures 2 and 3. As expected, the oxygen–oxygen pair correlation functions gOO(r) feature no significant differences for all of the order/disorder ice polymorph pairings. In addition, the oxygen–hydrogen pair correlation functions gOH(r) show only minute deviations, and only the hydrogen–hydrogen pair correlations gHH(r) differ appreciably. However, these differences are small, making it exceedingly difficult to distinguish between different variants of protonically ordered cubic ice based on the comparison of pair correlation functions determined experimentally from X-ray or neutron diffraction data.[36]

Figure 2

Partial radial distribution functions gOO(r) and gOH(r) for proton-disordered hexagonal and cubic ice as well as their proton-ordered counterparts at T = 170 K. Both the oxygen–oxygen as well as the oxygen–hydrogen radial distribution functions are nearly identical in all cases including proton-ordered and -disordered configurations. The curves labeled a, b, c, and d correspond to the proton orderings shown in Figure 1. Note that the peaks corresponding to the intramolecular OH distances lie outside of the range of the figure.

Figure 3

Partial radial distribution functions gHH(r) for proton-disordered hexagonal and cubic ice as well as their proton-ordered counterparts at T = 170 K. The labeling a, b, c, and d refers to the various proton-ordering patterns shown in Figure 1. Note that the peak corresponding to the intramolecular HH distances lies outside of the range of the figure.

More detailed structural information on proton ordering is encoded in vibrational spectra. For the calculations of the IR spectra using MD simulations, we have set up cubic and hexagonal lattices consisting of N = 1000 and 896 molecules, respectively. The time step was set to 2 fs, dipole–dipole time correlation functions were calculated for times up to 20 ps, and in total, each system was propagated for 40 ns using a Nosé–Hoover NVT integrator. Here, the Lennard-Jones cutoff and the Ewald summation real space cutoff were chosen to be 3σ. The densities for the different arrangements of hexagonal and cubic ice were set according to Table 1.

Figures 4 and 5 show the computed IR spectra for hexagonal and cubic ice, respectively. In general, the computed spectra are in reasonable agreement with experimental spectra. The broad librational band near 850 cm–1 is in good agreement with experimental spectra of both cubic and hexagonal ice (see the dashed line in Figure 3 of ref (56)). Also, in the far-IR (0–400 cm–1), the agreement is reasonable, even though there are some discrepancies especially regarding the band intensities. Hexagonal ice shows three broad bands in this area,[57] which are centered at 160, 230, and 370 cm–1. In our simulations, there are three broad bands centered at 70, 195, and 310 cm–1. The absorption spectra of the proton-disordered forms of hexagonal and cubic ice are clearly distinguishable from those of the various proton-ordered phases, which split up in several bands, the number and position of which depend on the particular proton-ordering pattern. Thus, these spectra provide a means of detecting the presence of proton order experimentally. IR spectra computed by density functional theory for proton-ordered cubic ice at T = 0 are shown in Figure 5 as vertical bars and confirm this conclusion. The length of each bar shown in Figure 5 corresponds to the dipole activity of the corresponding mode. Because the spectra were determined using a harmonic approximation, finite temperature broadening of the individual peaks is entirely absent. The overall structure of the ab initio spectra shares some important qualitative features with the spectra obtained from the MD simulations, although in the density functional theory calculations, the librational band (550–1100 cm–1) is shifted by about 100 cm–1 to higher frequencies. For instance, for structure (a), a single librational band at 930 cm–1 is predicted using density functional theory, whereas a 20 cm–1 broad peak at 856 cm–1 is observed in the MD simulation. Likewise, four peaks are predicted for structure (b), which correspond well to the four-peak structure observed in the MD simulations. A further noteworthy difference between the ab initio and the MD results is that the average frequency calculated ab initio shifts significantly to the blue with increasing ferroelectric ordering. We find average frequencies of 930, 900, 880, and 850 cm–1 for cubic ice with ordering (a) (ferroelectric), (b) (weakly ferroelectric), (c) (weakly ferroelectric), and (d) (antiferroelectric), respectively. In the MD simulations, the shift is smaller with frequencies of 849, 833, 826, and 816 cm–1 for the four phases. The difference is most likely related again to the neglect of polarization effects in the MD simulations using the TIP4P/ice model.

Figure 4

IR spectra obtained from MD simulations for disordered hexagonal ice at T = 170 K and two of its proton-ordered forms at T = 70 K. Vertical bars indicate the results of ab initio simulations at T = 0.

Figure 5

IR spectra for disordered cubic ice Ic and various of its proton-ordered forms (the labeling a, b, c, and d refers to the structures shown in Figure 1). Lines indicate results from the MD simulations carried out at T = 70 K, and vertical bars indicate the results of ab initio simulations at T = 0. For the disordered variant of cubic ice, also the spectrum obtained at T = 170 K from MD simulations is shown. The black dotted lines denote the experimental IR spectra of cubic ice at T = 160 K and after the cooling procedure at T = 70 K (same data as that shown in Figure 6).

Figure 6

Librational band of the IR spectra of cubic ice determined experimentally at T = 160 K (blue line) and after 80 h at T = 70 K (red line).

Experiments

As mentioned in section 1, ice Ih and ice Ic differ in their layer stacking but have an identical local structure. One would therefore expect that ice Ic shows a similar proton-ordering transition as ice Ih. The experimental protocol to proton order ice Ih involves bringing the sample to ∼50 K, which is thought to produce some ferroelectric ice XI seeds, and then waiting at ∼55–67 K for several days or weeks, which allows the ice XI domains to grow.[13−17] Typically, about 50–60% of the hexagonal ice sample orders using this procedure. Because the local ordering is identical in cubic and hexagonal ice, we anticipate that a similar protocol might allow for hydrogen ordering in cubic ice. To investigate this possibility, we have carried out experiments on cubic ice. As discussed in detail below, our observations indicate partial proton ordering detected by comparison of the measured IR spectra with the results of our simulations.

To carry out our experiments, we prepared an aqueous 0.1 M KOH solution containing 2 mol % D2O. The use of a small fraction of D2O allows for the observation of the decoupled OD-stretching mode in IR spectra, in addition to the coupled OH-stretching mode. KOH introduces substitutional point defects in the ice lattice, which increase the proton mobility and hence reorientational dynamics drastically. We used 0.1 M KOH in order to reach a saturation level of these substitutional point defects. In principle, also lower concentrations such as 0.01 M KOH might be suitable for reaching this enhancement in the dynamics. Cubic ice was prepared by spraying droplets of this solution of about 3 μm in diameter into a vacuum chamber containing a He-cryostated optical window, where the window was first kept at ∼77 K. This procedure is known to produce hyperquenched glassy water (HGW), which crystallizes to cubic ice upon heating to ∼160 K.[58,59] Cubic ice prepared in this way shows a comparably low number of hexagonal stacking faults and high cubicity index as judged from the intensity of the X-ray reflexes corresponding to the hexagonal faults.[59]

The crystallization from HGW to cubic ice was monitored by in situ FT-IR spectroscopy, and the librational band of cubic ice at ∼160 K is depicted as the blue curve in Figure 6. Cubic ice was then cooled to ∼60 K for a few hours and then heated to 70 K for about 80 h. The spectrum recorded after this procedure is depicted in red in Figure 6. It can clearly be noted that the half-width decreases for all observed bands. The librational band depicted in Figure 6 shows a decrease in fwhm from 251 to 196 cm–1, that is, to about 80%. In addition to the narrowing, also some structuring of the band is apparent such that additional Gaussians are required to fit the band. The band also loses intensity between 500 and 600 cm–1. From band decomposition analysis, two intense Gaussian components at 837 and approximately 905 cm–1, as well as a less intense band at 995 cm–1, are needed to explain the observed band shape. Note that a similar band narrowing was observed in IR spectra of the partial ordering of hydrogen atoms in ice Ih.[56] However, while there is no difference between the IR spectra of cubic ice and hexagonal ice at 160 K, shoulders emerge at different positions after partial hydrogen ordering, that is, the hydrogen order pattern differs between proton-ordered cubic ice and ice XI.

We now turn to the question whether the band changes seen in Figure 6 are a result of a simple thermal effect or whether proton ordering can be inferred from the data. Figure 7 shows the evolution of the librational band in KOH-doped ice Ic when changing the temperature step by step from 45 K to higher temperatures. The temperature broadening of the band is best illustrated in the magnified inset of Figure 7. It is quite small when increasing the temperature in the range 45–70 K (6 cm–1 per 25 K) but larger when changing the temperature from 70 to 75 K (4 cm–1 per 5 K). Above 75 K, the temperature broadening is again quite small (5 cm–1 per 25 K). This shows that an additional effect, besides thermal broadening, influences the change in the half-width of the librational band between 70 and 75 K, which we interpret to be proton disordering of partially ordered cubic ice. The occurrence of the proton-disordering temperature between 70 and 75 K in cubic ice seen in Figure 7 compares to the known proton-disordering temperature of 72 K in hexagonal ice.[10−17]

Figure 7

Libration mode of ice Ic doped with KOH at different temperatures, recorded after the following steps. Ice Ic doped with KOH was prepared by vitrification of the aqueous solution (200.2 g 0.1 M KOH(aq) + 4.9 g D2O) at 50 K using the technique of hyperquenching and subsequently crystallizing the amorphous deposit by bringing the sample holder for 0.5 h to 140 K. The following temperature protocol was then applied while continuously recording IR spectra: 14 h at 45 K; 2 h at 50 K; 20.5 h at 60 K; 7.5 h at 65 K; 22 h at 70 K; 23 h at 75 K; 7 h at 80 K; 2 h at 90 K; and 2 h at 100 K. After these steps, the temperature was decreased again to 70 K for 18.5 h. Selected spectra as marked are shown in the figure and magnified in the inset. The two thick lines show the librational mode after 22 h at 70 K (blue line), after the sample was kept at 45–65 K before, and after 18.5 h at 70 K (black line), after the sample was kept at 100 K before.

This implies that the phase boundary between ordered and disordered cubic ice is located at almost the same temperature as the phase boundary between ordered and disordered hexagonal ice. The two thick lines in Figure 7 compare two librational bands recorded for KOH-doped ice Ic, both determined at 70 K after having kept the sample several hours at 70 K: the narrower band (blue line) was obtained after several days of waiting at temperatures below 70 K and then heating to 70 K, whereas the broader band (black line) was obtained after cooling from 100 to 70 K, without any waiting time at temperatures below 70 K. In the case of the proton-ordering transition from hexagonal ice to ice XI, it was shown that the proton-ordered phase grows much faster at 70 K if it was nucleated before at temperatures below 70 K, for example, at 60 K,[10−17] whereas by cooling from 100 to 70 K, without a prior nucleation step at lower temperature, the proton ordering is much slower in case of KOH-doped hexagonal ice. We expect this to be similar in the case of KOH-doped cubic ice. The difference in half-width at half-maximum of about 6 cm–1 between the two spectra taken at 70 K seen in Figure 7 clearly demonstrates that this expectation is really seen in the spectra. Indeed, an effect other than motional narrowing is operative, which we explain as proton ordering when the proton-ordered phase was nucleated by keeping the sample for several days between 45 and 65 K. Note that also the peak intensities and baselines slightly shift from spectrum to spectrum, for example, when comparing the thick black and blue lines in Figure 7 or the black and red lines in Figure 8. The error bar on the changes in half-widths due to the ambiguities in half-height amount to at most ±0.5 cm–1. The observed changes in half-width of 6 and 10 cm–1 indicated in Figures 7 and 8, respectively, are much larger than the error bar of the method and, thus, are real effects associated with changes within the sample itself.

Figure 8

Libration mode of undoped ice Ic after very slow heating from 45 to 70 K (red line) compared with KOH-doped ice Ic after cooling from 100 to 70 K (black line, taken from Figure 7) and with KOH-doped ice Ic after very slow heating from 45 to 70 K (blue line, taken from Figure 7). All spectra were recorded at 70 K after several hours at 70 K.

We now turn to the question whether the proton ordering taking place at 70–75 K in cubic ice pertains to the cubic stacking sequences or to the hexagonal stacking faults. Depending on the route of preparation, cubic ice always contains more or less hexagonal stacking faults.[60−64] The hyperquenching technique of micrometer-sized droplets used here to prepare cubic ice[65] is advantegous because it is known to produce cubic ice of rather high cubicity (see Figure 1 in ref (59)) as compared to cubic ice prepared, for example, by vapor deposition or by a phase transition from high-pressure ices upon heating at ambient pressure. We emphasize that hexagonal ice Bragg peaks, for example, at 2θ = 34 or 40° (Cu–Ka), are not detected at all. The observed (noise level) intensity at 2θ = 34 or 40° is lower by about a factor of 1000 compared to that of the most intense (111) Bragg peak of cubic ice at 2θ = 24°. From this, we estimate that our cubic ice samples are contaminated with hexagonal ice levels of <0.1%. The level of hexagonal ice contamination as judged from the powder X-ray diffractograms remains <0.1% also for cubic ice samples prepared by hyperquenching a 0.1 M KOH solution rather than pure water.

Let us assume now that the additional narrowing of about 10 cm–1 observed in Figures 7 and 8 was caused by proton ordering of the two-dimensional layers of hexagonal stacking faults and that the dominating cubic stacking sequences were not proton-ordered. This would then imply that the band narrowing in a pure hexagonal ice crystal that only contains hexagonal layers would need to be >200 cm–1 upon its transition to ice XI. However, the band narrowing observed in pure hexagonal ice upon the ice Ih → ice XI transition at 72 K is also only about 10 cm–1.[56] In addition, one would expect that the band shape upon proton-ordering hexagonal stacking faults would resemble the band shape of ice XI. However, band decomposition shows that at <70 K, a band pattern appears that is clearly different from ice XI (see the discussion about Figure 6 above). That is, the experimental findings reported here rule out the possibility that the band narrowing is caused by proton ordering of hexagonal stacking faults in cubic ice. Instead, they suggest that the change in band shape is caused by the proton ordering of the cubic layers. It might be possible that both the cubic and the hexagonal stacking sequences may be proton-ordered. However, the number of hexagonal stacking faults is so small that we do not see any signs of proton ordering to ice XI in our spectra, that is, their influence on the IR band shape is negligible.

Next, we consider the question whether proton ordering of cubic ice is also possible in the case of undoped cubic ice samples. Figure 8 shows a comparison of the two librational bands recorded at 70 K on KOH-doped cubic ice (black and blue traces, both taken from Figure 7) and on undoped cubic ice (red trace), which was kept several days at 45–65 K before. Clearly, the librational band is narrower in the case of KOH-doped ice. In this case, the difference in half-width is even about 10 cm–1. This shows that KOH doping is required to achieve proton ordering in cubic ice, that a smaller degree of proton ordering is also reached when cooling directly from 100 to 70 K, without the nucleation step at 45–65 K, whereas in the case of undoped cubic ice, proton ordering cannot be inferred from the data. Further systematic studies about the dependence of the half-width of the librational band on the thermal history using much longer waiting times (on the order of months) are necessary to find conditions that might allow for a higher degree of proton ordering or even to access the fully proton-ordered state of cubic ice. Also, an investigation of the influence of other dopants (e.g., NH3, HCl, etc.) is of interest in this context.

Discussion

Assuming that the band structuring and narrowing that is observed in addition to the pure thermal narrowing in our spectra indeed arises from partial hydrogen ordering in cubic ice, we may now compare the observed band with the predicted spectra of Figure 5. For protonically disordered cubic ice, the predicted librational band centered at ∼800 cm–1 coincides almost perfectly with the experimentally observed librational band both in position and width, which is about fwhm ≈ 200 cm–1. Our MD simulations predict that depending on the particular type of proton ordering, this band splits into up to six bands, all of which show a fwhm of about ∼20 cm–1. That is, the fwhm of the ordered forms is about 10% of the fwhm in the case of the disordered form. In the experiment, the fwhm in the possibly ordered form of cubic ice is reduced to about 80%, thus suggesting that the degree of hydrogen ordering achieved experimentally is no more than 20%.

Band deconvolution suggests that three bands contribute to the experimental band shape. The main peak at 837 cm–1 is most likely related to the predominant disordered cubic ice and shifts only little compared to the original disordered cubic ice phase (816 cm–1). This small shift of 20 cm–1 might be related to the lower temperature of 70 K compared to 160 K. While usually there is a blue shift when lowering temperature, here we observe a red shift, which is probably related to the anomalous temperature dependence of ice[1] and may also result in positive Grüneisen parameters for several modes.[66] Furthermore, in the experiment, a distinct shoulder at 905 cm–1 is observed accounting for about one-third of the integrated intensity.

This band structure of the partially ordered phase seems to be consistent with a more pronounced ferroelectric order, which both in the DFT calculations as well as in the MD simulations leads to a pronounced blue shift compared to antiferroelectric order. In particular, for the ferroelectric structure (a), DFT predicts a strong single peak at 930 cm–1, while the vibrational band of the antiferroelectric order (d) is shifted to the red by about 80 cm–1 on average. In the TIP4P/ice simulations, the dominant peak for cubic ice (a) as well as (b) is located at a lower frequency of 856 cm–1, a difference of 70 cm–1.

The experiment also shows a weak but clearly resolved shoulder at 995 cm–1, approximately 90 cm–1 above the main peak arising as a result of the conjectured partial proton order. This side peak is not accounted for by the ferroelectrically ordered structure (a), which shows only one single peak. Such side peaks are, however, present in both the cubic ordered variants (b) and (c), with structure (b) showing very similar ratios between the main peak and the higher-frequency shoulder as in the experiment. In fact, the magnitude of the two band splittings is about 70 and 90 cm–1 in the experiment and thus quite similar to the splittings of 90 and 90 cm–1 predicted for the ordered structure (b) (see Figure 5). The other three ordered variants of cubic ice show quite different patterns and splittings. This might suggest that the transformation occurs into structural variant (b) and not the lowest-energy structure ordered ice (a), so that the experimentally observed ordered form does not necessarily need to correspond to the thermodynamically most stable form. This is plausible because kinetics is known to play an important role in the proton-ordering transitions of ice, as, for instance, observed in ice XIV.[67] It is also possible that the experimental spectrum might develop into the spectrum of ordered ice (a) with increasing time and increasing proton ordering instead of developing the peak splittings predicted for structure (b).

Interestingly, the ordered variants of cubic ice also show shifts and/or splittings of the acoustic and optical modes at <400 cm–1 compared to disordered cubic ice according to our simulations (see Figure 5). That is, investigation of acoustic and optic modes in the future by means of IR or Raman seems to be a promising tool for critically testing the type of proton ordering obtained experimentally in cubic ice.

The temperature dependence of the librational band shape and width suggest that the proton order–disorder transition in cubic ice takes place at about 70–75 K, which is very similar to the proton order–disorder temperature of 72 K observed in the case of hexagonal ice. Furthermore, the degree of proton ordering and narrowing of the bands at 70 K is clearly enhanced when employing a “nucleation” step at temperatures of 45–65 K. Whereas we observe in our experiments on KOH-doped cubic ice at T < 75 K a band narrowing in addition to the narrowing caused by motional narrowing, we do not observe such an additional narrowing in the case of undoped, pure cubic ice on the time scale of days.

In summary, we have studied the proton-ordering transition of ice Ic using a combination of FT-IR spectroscopy experiments and computer simulations. We find that subjecting a sample of ice Ic to an experimental protocol similar to that used before to proton order hexagonal ice causes significant changes in the librational band both in shape and position. A comparison of the experimental IR spectra with theoretical spectra obtained from computer simulations indicates that these changes are likely due to partial proton ordering in the ice Ic sample occurring in the temperature range of 70–75 K. On the basis of a comparison of the computed and measured librational spectra, no unique assignment of the type of protonic order is currently possible. While the intensity loss at low frequencies suggests a partial ordering into a ferroelectric structure, which has the lowest energy, the number and relative positions of the peaks obtained from deconvolution of the experimental data point to a weakly ferroelectric structure with a slightly higher energy. Further experiments and simulations will be necessary to resolve this issue.