Classical and quantum theories of proton disorder in hexagonal water ice

It has been known since the pioneering work of Bernal, Fowler, and Pauling that common, hexagonal (Ih) water ice is the archetype of a frustrated material: a proton-bonded network in which protons satisfy strong local constraints (the “ice rules”) but do not order. While this proton disorder is well established, there is now a growing body of evidence that quantum effects may also have a role to play in the physics of ice at low temperatures. In this paper, we use a combination of numerical and analytic techniques to explore the nature of proton correlations in both classical and quantum models of ice Ih. In the case of classical ice Ih, we find that the ice rules have two, distinct, consequences for scattering experiments: singular “pinch points,” reflecting a zero-divergence condition on the uniform polarization of the crystal, and broad, asymmetric features, coming from its staggered polarization. In the case of the quantum model, we find that the collective quantum tunneling of groups of protons can convert states obeying the ice rules into a quantum liquid, whose excitations are birefringent, emergent photons. We make explicit predictions for scattering experiments on both classical and quantum ice Ih, and show how the quantum theory can explain the “wings” of incoherent inelastic scattering observed in recent neutron scattering experiments [Bove et al., Phys. Rev. Lett. 103, 165901 (2009)]. These results raise the intriguing possibility that the protons in ice Ih could form a quantum liquid at low temperatures, in which protons are not merely disordered, but continually fluctuate between different configurations obeying the ice rules.

Figure 1

Crystal structure of hexagonal (Ih) water ice. Water ice can be viewed as a frozen configuration of water molecules, satisfying the Bernal-Fowler “ice rules” [1, 2], in which each oxygen (red sphere) forms two short, covalent bonds, and two long, hydrogen bonds with neighboring protons (white spheres). Oxygen atoms form an ordered lattice, belonging to the hexagonal space group $P{6}_3/mmc$, with a four-site primitive unit cell. Protons do not show any long-range order. (a) Structure viewed perpendicular to the hexagonal symmetry axis (the crystallographic $c$ axis). (b) Structure viewed along the hexagonal symmetry axis.

Figure 2

Collective quantum tunneling between different proton configurations satisfying the ice rules in hexagonal (Ih) water ice. The crystal structure of ice Ih contains two distinct types of hexagonal plaquette, each containing six oxygen atoms (red spheres). Tunneling between different ice configurations is mediated the by correlated hopping of protons (white spheres) around such a plaquette. Where protons form an alternating sequence of short (S) and long (L) bonds (e.g., S-L-S-L...-S-L), it is possible to tunnel to another, degenerate, ice configuration in which long and short bonds are interchanged (i.e., L-S-L-S...-L-S) [18, 19, 20]. (a) Tunneling on hexagonal plaquette of “type I” in the plane perpendicular to the hexagonal symmetry axis. This process has a matrix element $g_1$ in our model for proton dynamics ${\mathcal{H}}_{\mathsf{\text{tunneling}}}^{\mathsf{\text{hexagonal}}}$ [Eq. (4)]. (b) Tunneling on hexagonal plaquette of “type II,” connecting different hexagonal planes. This process has a matrix element $g_2$.

Figure 3

Different types of bond within the ordered oxygen lattice of hexagonal (Ih) water ice. The Ih lattice can be viewed as a stack of buckled, honeycomb lattices, composed of bonds which are symmetric under inversion about the center of the bond. These honeycomb layers are connected by bonds running parallel to the hexagonal symmetry axis, which are symmetric under reflection in the plane perpendicular to the bond. This lattice is bipartite, i.e., it may be divided into two sublattices, here colored red and blue. (a) Center-symmetric bond, viewed from a direction perpendicular to the bond. (b) Center-symmetric bond, viewed along the bond direction. (c) Mirror-symmetric bond, viewed from a direction perpendicular to the bond. (d) Mirror-symmetric bond, viewed along the bond direction.

Figure 4

Flux representation of proton configurations in hexagonal (Ih) water ice. (a) Proton configuration, including a hexagonal plaquette of type II (cf. Fig. 2). (b) Equivalent proton configuration represented using arrows. The displacement of protons (white spheres) from the midpoint of each oxygen-oxygen bond can be mapped to an arrow on that bond. The ice rules require that two “in” arrows and two “out” arrows meet at each vertex of the lattice (oxygen atom). Where arrows form a closed loop, it is possible to tunnel between different proton configurations satisfying the ice rules, by reversing the sense of all arrows on that loop. The letters $A,B,C,D$, and color coding, indicate the convention for labeling oxygen sublattices adopted in Sec. 3.

Figure 5

Correlation of long and short proton bonds in a classical model of hexagonal (Ih) water ice. The structure factor $S_{\text{Ising}}\left(\mathbf{q}\right)$ [Eq. (29)], for the Ising variable $\sigma$ [Eqs. (7)–(10)], is plotted in three orthogonal planes in reciprocal space. Near to zone centers, correlations are well described by a combination of pinch-point singularities, reflecting the algebraic correlations of the field ${\mathbf{P}}_{+}\left(\mathbf{q}\right)$ [cf. Eq. (21)], and smooth features, reflecting the short-range correlations of the field ${\mathbf{P}}_{-}\left(\mathbf{q}\right)$ [cf. Eqs. (22) and (23)], as discussed in Sec. 3a. Calculations were performed using the method outlined in Appendix pp3, in which the ice-rule constraints are written as orthogonality conditions in Fourier space [57]. Reciprocal-lattice vectors are indexed to the orthorhombic unit cell defined in Appendix pp1, following the conventions of Nield and Whitworth [52].

Figure 6

Prediction for diffuse scattering of neutrons or x rays from protons in a classical water ice described by the Bernal-Fowler “ice rules.” The structure factor $S_{\mathsf{\text{proton}}}^{\mathsf{\text{diffuse}}}\left(\mathbf{q}\right)$ [Eq. (45)] is plotted in three orthogonal planes in reciprocal space. The ice rules manifest themselves as both “pinch points” [singular features in scattering visible at, e.g., ${\mathbf{Q}}_{\mathsf{p},{\mathsf{H}}^{+}}^{*}=(0,0,4)$ in (c)] and broad, asymmetric features at zone centers visible at, e.g., ${\mathbf{Q}}_{\mathsf{m},{\mathsf{H}}^{+}}^{*}=(2,0,3)$ in (c). Results are shown for the theory described in Sec. 3c and Appendix pp3, with reciprocal-lattice vectors indexed to the orthorhombic unit cell defined in Appendix pp1, following the conventions of Nield and Whitworth [52].

Figure 7

Prediction for inelastic scattering from disordered protons in a quantum model of ice Ih at zero temperature ($T=0$). Results for the dynamical structure factor $S_{\mathsf{coh}}^{{\mathsf{H}}^{+}}(\mathbf{q},\omega )$ [Eq. (E24)] are shown in three, orthogonal, planes in reciprocal space. The “photons” of the lattice-gauge theory ${\mathcal{H}}_{\mathsf{\text{U}}\left(\mathsf{1}\right)}$ [Eq. (66)] are visible as gapless, linearly dispersing excitations at long wavelength. These photons are birefringent, with a dispersion which depends on the polarization of the photon, except along the optical axis $\mathbf{z}$. As a result, two bands of photons are visible on the path $\mathrm{\Gamma }-M-K-\mathrm{\Gamma }$, while only one appears on the path $\mathrm{\Gamma }-A-L-H-A$. Additional spectral weight at higher energies is associated with gapped, exponentially correlated excitations of the gauge field. The dynamical structure factor for coherent scattering from protons $S_{\mathsf{coh}}^{{\mathsf{H}}^{+}}(\mathbf{q},\omega )$ [Eq. (E24)] was calculated using the $\mathsf{\text{U}}\left(\mathsf{1}\right)$ lattice-gauge theory ${\mathcal{H}}_{\mathsf{\text{U}}\left(\mathsf{1}\right)}$ [Eq. (66)]. Results are shown for the parameter set given in Eq. (68), and were convoluted with a Gaussian of $\text{FWHM}=0.035\sqrt{\mathcal{U}\mathcal{K}}$ to mimic the effect of experimental resolution. The color scale shows the intensity of the modes as they would appear in an inelastic scattering experiment, in the Brillouin zone centered on $(h,k,l)=(0,2,4)$.

Figure 8

Prediction for coherent, quasielastic scattering from protons in a quantum model of ice Ih at zero temperature ($T=0$). The pinch points associated with the ice rules (cf. Fig. 6) are eliminated by quantum fluctuations (cf. Ref. [39]). At finite temperatures, these pinch points will be restored with a weight linear in $T$ [37]. Results are shown for the energy-integrated, equal-time structure factor $S_{\mathsf{coh}}^{{\mathsf{H}}^{+}}(\mathbf{q},t=0)$ [Eq. (79)], calculated within the lattice-gauge theory ${\mathcal{H}}_{\mathsf{\text{U}}\left(\mathsf{1}\right)}$ [Eq. (66)], for the same parameter set [Eq. (68)] as Fig. 7. Reciprocal-lattice vectors are indexed to the orthorhombic unit cell defined in Appendix pp1, following the conventions of Nield and Whitworth [52].

Figure 9

Prediction for inelastic, incoherent scattering from protons in quantum ice Ih. (a) Prediction of the lattice-gauge theory ${\mathcal{H}}_{\text{U}\left(1\right)}$ [Eq. (66)], at a temperature of $T=5\text{K}$, for the same parameters set used in discussing coherent inelastic scattering (cf. Figs. 7 and 8). (b) Prediction of the lattice-gauge theory, convoluted with a Gaussian of FHWM $0.07\text{meV}$ to represent finite experimental resolution. (c) Prediction of the lattice-gauge theory (blue), combined with the elastic contribution to the incoherent scattering (yellow). Both have been convoluted with a Gaussian to represent finite experimental resolution. The combined line shape shows inelastic “wings,” similar to those observed in experiments on ice Ih by Bove et al. [18]. Parameters are given in Eqs. (68) and (86), with details of calculations in Appendix pp7.

Figure 10

Labeling convention for the eight sets of oxygen-oxygen bonds and eight sets of six-link plaquettes unrelated by translational symmetry. This is the convention employed in diagonalizing the lattice-gauge theory described in Sec. 4 and Appendix pp5. (a) Color coded such that bonds related by translational symmetry are the same color. The bonds labeled 0 and 4 are mirror-symmetric bonds aligned along the $\mathbf{z}$ axis and the remaining bonds are center symmetric (cf. Fig. 3). (b) Shows plaquettes of type I with matrix element $g_1$. (c)–(e) Shows plaquettes of type II with matrix element $g_2$.

Figure 11

Identification of the excitations of the lattice-gauge theory ${\mathcal{H}}_{\mathsf{\text{U}}\left(\mathsf{1}\right)}$ [Eq. (E1)], with quantized fluctuations of the fields ${\mathbf{P}}_{+}$ and ${\mathbf{P}}_{-}$. (a) Dynamical structure factor of the uniform polarization ${\mathbf{P}}_{+}$ [Eq. (E25)], calculated within the lattice-gauge theory. (b) Equivalent results for the dynamical structure factor of the staggered polarization ${\mathbf{P}}_{-}$. For $q\to 0$, the correlations of ${\mathbf{P}}_{+}$ are directly associated with the gapless, birefringent photons of the lattice-gauge theory, while the correlations of ${\mathbf{P}}_{-}$ are associated with its gapped optical modes. This identification breaks down at shorter wavelengths, where the correlations of ${\mathbf{P}}_{+}$ and ${\mathbf{P}}_{-}$ have contributions from all modes of the lattice-gauge theory. Calculations were carried out for the symmetric choice of parameters $\mathcal{U}={\mathcal{U}}^{\prime }$ and $\mathcal{K}={\mathcal{K}}^{\prime }$.

Figure 12

Ising structure factor $S_{\mathsf{\text{Ising}}}\left(\mathbf{q}\right)$ [Eq. (29)] for a classical model of ice Ih, calculated from classical Monte Carlo simulation [(a)–(c)] and from lattice theory [(d)–(f)]. The theory calculation is performed using the method outlined in Appendix pp3, which is based on the method described for spin ice in Ref. [57]. This method consists in writing the ice-rule constraints as orthogonality conditions in Fourier space. In the long-wavelength limit, these conditions reduce to those obtained from the continuum field theory presented in Sec. 3a [Eqs. (17) and (18)]. Near the zone centers, the correlations are well described by a combination of pinch-point singularities, reflecting the algebraic correlations of ${\mathbf{P}}_{+}\left(\mathbf{q}\right)$ and smooth features reflecting the short-ranged correlations of ${\mathbf{P}}_{-}\left(\mathbf{q}\right)$.

Figure 13

Comparison between the predictions of the lattice-gauge theory ${\mathcal{H}}_{\mathsf{\text{U}}\left(\mathsf{1}\right)}$ [Eq. (66)], and variational Monte Carlo simulation of ${\mathcal{H}}_{\mathsf{\text{tunneling}}}^{\mathsf{\text{hexagonal}}}$ [Eq. (4)] for the “Ising” structure factor $S_{\mathsf{\text{Ising}}}(\mathbf{q},t=0)$ [Eq. (F2)], for the parameter sets given in Eqs. (F3) and (F4). There is excellent agreement between the two methods, validating our description of the proton liquid in water ice Ih in terms of a quantum $\mathsf{\text{U}}\left(\mathsf{1}\right)$ lattice-gauge theory.