First-principles momentum distributions and vibrationally corrected permittivities of hexagonal and cubic ice

Three-dimensionally-resolved proton momentum distributions and end-to-end distributions have been calculated for hexagonal and cubic water ice. First-principles quantum nuclear wave functions have been used to investigate the impact of vibrational anisotropy, anharmonicity, proton-and stacking-disorder, temperature, and pressure on these distributions. Moreover, the eﬀects of vibrations on the electronic density in hexagonal ice are shown to lead to a 5 % vibrational correction with respect to the static-lattice optical permittivity, and proton-disorder is found to be crucial in explaining its experimentally observed temperature dependence.

Three-dimensionally-resolved proton momentum distributions and end-to-end distributions have been calculated for hexagonal and cubic water ice. First-principles quantum nuclear wave functions have been used to investigate the impact of vibrational anisotropy, anharmonicity, protonand stacking-disorder, temperature, and pressure on these distributions. Moreover, the effects of vibrations on the electronic density in hexagonal ice are shown to lead to a 5 % vibrational correction with respect to the static-lattice optical permittivity, and proton-disorder is found to be crucial in explaining its experimentally observed temperature dependence.
Hydrogen-bonded materials play an important role across many fields of science. The importance of water ice, in particular, spans disciplines ranging from astrophysics to biology. Properties of hydrogen-bonded materials are linked to those of the hydrogen bond [1,2], and are strongly affected by the quantum nuclear motion of the light protons [3][4][5][6][7]. Anharmonic quantum nuclear effects (QNE) play a key role in ice. For example, they stabilise hexagonal ice (Ih) with respect to cubic ice (Ic), which ultimately leads to the hexagonal shape of ice crystals [7]. They also play an important role in the anomalous thermal expansion of ice Ih [8], proton/deuteron isotopic effects [9][10][11], and in shifts in infra-red and other vibrational spectra [7]. Quantum zero-point (ZP) and thermal motion of the protons has a large impact on the electronic properties of ice [12,13]. Here we study the effects of vibrational motion on the electronic density and permittivity of ice using first-principles densityfunctional-theory (DFT) methods. Quantum ZP motion dominates the equilibrium proton dynamics of ice up to the melting temperature [12]. Proton dynamics therefore provide a rare direct probe of their quantum nature. We have compared the results of our first-principles calculations of position and momentum distributions with experimental data, which provides a stringent test of computational descriptions of ice and the hydrogen bond.

I. PROTON REAL SPACE AND MOMENTUM DISTRIBUTIONS
The proton radial distribution function (RDF) provides a real-space measure of proton dynamics, whereas the momentum distribution function (MDF) gives complementary reciprocal-space information [14]. Experimentally, the RDF can be accessed by inversion of the neutron [15] and x-ray structure factors [16], while neutron Compton scattering and deep inelastic neutron scattering (DINS) [17,18] can probe the MDF. Path-integral molecular dynamics (PIMD) simulations have been the method of choice for calculating both the RDF [19][20][21] and MDF [22][23][24] and have provided important insights into the RDF and MDF signatures of the interactions of the protons with their environment. However, they have also demonstrated that some of the most commonly used empirical models of water do not accurately reproduce the experimental MDF. Experimental data is currently limited to the spherically averaged MDF and is insufficient to disentangle the effects of environmental anisotropy and vibrational anharmonicity [25]. This gives rise to severe difficulties in interpreting experiments on nano-confined and supercooled water [26,27]. Such experiments highlight the need for accurate ab initio simulations of proton distribution functions, such as the path-integral Car-Parinello molecular dynamics (PICPMD) studies of Lin et al. [25] and Flammini et al. [18].
PIMD studies, including those of the PICPMD flavour have their limitations: (1) they are expensive, which limits the number (and sizes) of systems, temperatures, and pressures that can be studied. (2) The effects of anharmonicity cannot be extracted directly, as PIMD inherently simulates the behaviour of the anharmonic system. (3) At its core PIMD is a sophisticated phase space sampling method in which statistical uncertainty limits the accuracy with which (complicated) quantities, such as the spatially-resolved MDF can be determined. To overcome these limitations, we calculate the MDF directly from the fully-anharmonic nuclear wave function that is obtained using the vibrational self-consistent field (VSCF) approach of Monserrat et al. [28]. This provides the unprecedented level of microscopic insight required to disentangle the effects of anisotropy, vibrational anharmonicity, pressure, temperature, and stacking-and proton-disorder. Our analytic description of the vibrational wave function also greatly facilitates calculations of accurate three-dimensionally-resolved MDFs. Moreover, the VSCF approach naturally lends itself to the use of non-diagonal supercells [29], which enable extensive sampling of the vibrational Brillouin Zone (BZ) without the need for large supercell simulations. This renders the VSCF approach far less expensive than PIMD methods when points in the vibrational BZ other than the Γ-point need to be sampled for converged results. Notably, calculations using 64-molecule simulation cells show that Γ-point calculations using eight-molecule simulation cells suffice in the context of proton RDFs and MDFs in ice Ih and Ic.
Computational methods. We use DFT as implemented in the Castep plane-wave pseudo-potential package [30] (version 7.03) with the PBE exchangecorrelation (xc) functional [31], which reproduces the experimental lattice parameters of ice Ih and Ic rather well [7]. Geometry optimisations and the mapping of the Born-Oppenheimer potential energy surface (PES) were performed with a plane-wave energy cut-off of 1600 eV, a Monkhorst-Pack k-point grid [32] spacing of less than 2π × 0.04Å −1 , and on-the-fly generated ultrasoft pseudopotentials [33]. The resulting energy differences between frozen-phonon configurations are converged to below 10 −4 eV/H 2 O, the atomic positions are converged to within 10 −5Å , and the residual forces to within 10 −4 eV/Å. For an N -atom system, the 3N harmonic vibrational normal modes (n, K) with branch index n, first vibrational BZ momentum K, and frequencies ω nK are calculated using a finite displacement approach [34]. The harmonic approximation is used to define normal mode coordinates q nK , which measure the collective atomic displacements along the normal modes (n,K). They are related to the Cartesian displacement of the nucleus i with mass M i from its equilibrium position, u i , by where w −KnIµ is the displacement pattern of mode (n, K). The momentum K is suppressed in the following as we restrict ourselves to Γ-point calculations.
The light mass of the protons, M p = 1837.3621238166 au, leads to large displacements due to quantum zero-point motion, implying that anharmonicity in the PES is important. In ice Ih and Ic anharmonic vibrations can be calculated within the independent-mode approximation due to the small coupling between normal modes [7]. Accordingly, the PES is expressed as a function of q = (q i , . . . , q 3N ), expanded in independent-mode terms V (1) (q n ), pairwise couplings V (2) (q n , q n ′ ), etc., and then truncated after (2) The 3N anharmonic independent-mode potentials V (1) (q n ) are mapped up to large amplitudes of four times the harmonic root-mean-square (RMS) displacements and fitted using cubic splines. Mapping each V (1) (q n ) using 11 equally spaced points was found to lead FIG. 1. Comparison of the spherically-averaged momentum distribution np(p) from experiment (yellow [23], orange [17], and red [18]) and path-integral Car-Parinello molecular dynamics (PICPMD) simulations (blue) [18] to VSCF vibrational calculations (black). The inset compares the mean force f (x) with error bars inferred from neutron Compton scattering (red) [18] to PICPMD data (blue) [18] and VSCF calculations. The deviation of f (x) from linearity indicates that the proton potential exhibits substantial anharmonicity. The uncertainties in the experimental distributions arise from fitting models to data with finite experimental resolution. The PICPMD data exhibits small statistical errors, while the VSCF results exhibit negligibly small errors due to the discretisation of n(r) and n(p).
to converged results. The vibronic Schrödinger equation is solved following the vibrational self-consistent field method described in Ref. [28]. The anharmonic vibrational wave function is expanded as a product, of single-particle anharmonic eigenstates, |ψ Sn (q n ) , each using a basis of 25 simple harmonic oscillator eigenstates. S denotes the vibrational eigenstate whose elements S n label the states of the independent modes. Notably, the choice of xc-functional affects nuclear vibrations in ice Ih and Ic very little for fixed unit cell parameters [7].
Comparison with experiment and PICPMD results. Using the anharmonic (unless explicitly stated) vibrational wave function, the distribution n i (p) of the momentum p of proton i can be calculated as the Fourier transform of the one-body reduced density ma-  (7) to np(p) at 268 K, while Tp denotes direct, non-parametric estimates of the kinetic energy of the protons. σ and σ ⊥ measure the widths of n(p) parallel and perpendicular to the O-H direction, averaged over proton-orderings. Errors arise from fitting Eq. 8 to noisy data due to the experimental resolution (Exp), statistical uncertainty from finite sampling (PICPMD), and the discretisation of n(r) and n(p) (VSCF), respectively.
where (d 3 r j =i ) denotes integration over the atomic positions r j of all 3N atoms in the simulation cell except atom i. In the following we refer to as the "end-to-end distribution" in analogy with openended path-integral simulations [35], and the absence of the index i denotes quantities that have been averaged over all protons. In practice n i (r) and n i (p) are both descretised in real and reciprocal space, respectively. Measurements of the spatially-resolved MDF, n(p), are not yet available for ice, but the spherically averaged distribution n p (p) ≡ dθ dφ p 2 n(p) has been measured, for example, by Reiter et al. [17] and Senesi et al. [36]. Their data for n p (p), the anisotropy of n(p), and the kinetic energy of the protons provide a measure of the accuracy of computational approaches. The spherically averaged MDF, n p (p), is conventionally described by where σ and the dimensionless a n are fitting parameters and the H n are the usual Hermite polynomials. The widths σ obtained in this work and listed in Table I agree well with recent experimental [18] and PICPMD [25] data. We find that the kinetic energy of the protons T p , which can be calculated as T p = ∞ 0 p 2 /(2M )n p (p)dp or (within the harmonic approximation) via T p = 3 2 /(2M )σ 2 , match the experimental and PICPMD data [18,25].
The insensitivity of n p (p) to anisotropy in the momentum distribution popularised the mean force f (x) as a more sensitive probe of the local environment of a proton [35]. Experimentally, f (x) can be extracted via its relationship with the Fourier transform of the neutron Compton profile. Its slope at x = 0 determines the curvature of the proton potential around the equilibrium position and thereby provides a second measure of the kinetic energy of the protons. A comparison of the mean force obtained using the VSCF vibrational method to that from PICPMD simulations [18] is shown in the inset of Fig. 1.
The anisotropy of n(p) is conventionally described in terms of the widths of the three-dimensional anisotropic Gaussian fit to n i (p), where d and the widths σ α are fitting parameters. The O-H· · · O axis is characterised by σ , while the two normal axes are typically characterised by a single parameter σ ⊥ . The differences between the in-and out-of-molecular plane environments of each proton imply that the two normal axes are inequivalent, and the widths of n i (p) along the normal axes indeed differ by about 10%. Disregarding this, Table I demonstrates that the MDF from the VSCF method is about as good as the PICPMD counterparts, provided vibrational anharmonicity is accounted for. It also demonstrates good agreement with the more recent NCS data from Ref. [18], with the exception of the anisotropy of n(p). The differences in σ and σ ⊥ are likely due to the finite number of proton-orderings considered.
Measurements of the spatially-resolved momentum distribution in KH 2 PO 4 [37], which in its crystalline form is used in optical modulators and for non-linear optics such as second-harmonic generation, and Rb 3 H(SO 4 ) 2 [38] suggest that n i (p) in Ih should become accessible in the future with improved experimental techniques.
End-to-end distributions. In line with Eq. (8), the end-to-end distribution, as defined in Eq. (6), of a proton i in proton-disordered ice is well approximated by an anisotropic Gaussian, n i (r) ≈ exp(−1/2r T C −1 i r) [25], where C i,α,β = r i,α r i,β is the correlation matrix of proton i. The principal axes of the anisotropic Gaussian form reflect the molecular orientation of the particular proton.
However, in proton-ordered ice the local environments of the protons are generally inequivalent, which is reflected in the end-to-end distribution of individual protons. To highlight the anisotropy of n i (r) and the role of the local proton environment, we factorise n i (r) into an isotropic freeparticle contribution and an anisotropic environmental componentñ i (r), n i (r) ≡ñ i (r) exp(M p k B T r 2 ), where M p is the proton mass and T is the temperature. Fig. 2 (a) shows thatñ i (r) for an individual proton in Ih Cmc2 1 -the lowest energy proton-ordering known in Ih -deviates substantially from the purely anisotropic Gaussian form and reflects, in particular, the positions of the next-nearest-neighbour protons. Moreover, in a particular proton-ordered configuration of Ih, such as Cmc2 1 , not all possible environments are realised, so that even the proton-averagedñ(r) shown in Fig. 2 (b) deviates from the hexagonal symmetry expected from a anisotropic Gaussian form [25]. Averaging over the 16 eight-molecule proton-orderings of Ih used by Hirsch and Ojamae [39] as shown in Fig. 2 (c) recovers the hexagonal symmetry expected from proton-disordered Ih.
Stacking-and proton-disorder and vibrational anharmonicity. As shown in Fig. 3 (a) the proton momentum distributions of different proton-orderings of Ih are indistinguishable. The differences between the structures of Ih and Ic are very small but are reflected even in the spherically averaged MDFs. The stronger localisation of the protons in Ic evidenced by 2 % smaller root-meansquare (RMS) displacements [7] manifests itself in a MDF that is around 2 % wider in Ic than in Ih. Real ice "Ic" samples are typically stacking-disordered [40,41]. They contain both Ih and Ic, which differ only in the stacking of molecular layers and which form the end members of the infinite set of stacking-disordered structures generated by introducing stacking faults into the ground-state structure Ih. The MDFs of Ih and Ic thus limit the effects of stacking-disorder on the spherically averaged proton MDFs of real ice samples.
The RMS proton displacements of Ref. [7], moreover, indicate that vibrational anharmonicity localises the protons in Ih by around 2 % with respect to the harmonic proton density distribution, suggesting that vibrational anharmonicity should widen the distribution n(p). Overall, the differences between the MDFs of different proton-orderings of Ih and Ic are too small to distinguish in experiment, but the effects of vibrational anharmonicity, while still small, are already of comparable size to the uncertainty in recent NCS data [18]. This suggests that it should be possible to discern them in future experiments.
Temperature and pressure dependence. The effects of vibrational anharmonicity and stacking-disorder on n p (p) are small, begging the question of the roles played by temperature, thermal expansion and pressure.
In ice only the softest modes, i.e., pseudo-translations, are thermally excited up to melting. Their contribution to proton motion is small, so that any temperature dependence must arise from thermal expansion. The volume expansion of Ih between zero temperature and melting is around 2 %. However, neither n p (p) nor n(p) exhibit any discernible temperature dependence across this temperature range.
Conversely, external pressures of up to ∼ 200 MPa lead to volume compressions of up to 1.8 %, at which point the Ih-II or Ih-III transition occurs. However, n p (p) or n(p) show no significant pressure dependence, even though the RMS proton displacements decrease by up to 0.5 %.
The insensitivity of n p (p) to pressure and temperature arises from the same principle: while pseudo-translations and librational modes soften upon expansion, the molecular modes harden, leading to cancellation of their effects on proton motion. Evidence of this is provided in the supplementary information (SI) [42].

II. VIBRATIONAL RENORMALISATION OF ELECTRONIC PROPERTIES
Studies of vibrational corrections to band gaps of bulk ice Ih and Ic [12] and their surface band gaps [13] demonstrate that the electronic density is substantially renormalised by quantum nuclear motion. We make use of the explicit form of the anharmonic nuclear wave function provided by the VSCF approach and the adiabatic approximation to sample the electronic density, Born effective charges, and permittivity using Monte Carlo methods.
Within the Born-Oppenheimer approximation, expectation values of vibrationally renormalised properties, O(T ), at temperature T may be written as is the partition function,Ô g (q) is the value of observable O for a frozen-phonon structure with atomic positions q = (. . . , q nk , . . .). The summation over vibrational eigenstates S includes the ground-state and excited states.
Here we use Monte Carlo sampling of O(T ) with N s frozen-phonon structures, q i , randomly drawn from the vibrational density, as 1/Z(T ) S Ψ S (q) 2 e −βE S , which gives We initially sample the electronic density, Born effective charges, and permittivity using the harmonic (har) vibrational density at a high temperature T ′ = 260K, and where β ′ = 1/(k B T ′ ). This reweighting approach is accurate, because the target vibrational density distributions are narrower than the harmonic distribution at 260 K, from which the frozen-phonon samples are drawn originally.
In practice at least 350 frozen-phonon configurations q i were sampled to determine various properties for each proton-ordering, reducing statistical uncertainties to less than 0.5%. By considering the same frozen-phonon configurations for different functionals, five configurations using the hybrid HSE06 functional [43,44] and about 20 configurations for each additional semi-local functional other than PBE are sufficient to gain insight into the role of the xc-functional. More detail regarding the choice of functional can be found in the SI. Delocalisation of the valence electron density due to quantum nuclear motion has a striking effect not only on the valence band maximum in ice and its electronic quasiparticle band gap [12], but also on the electronic and nuclear polarisability of ice. The nuclear polarisability, in particular, is affected by the renormalisation of the Born effective charges, Z ⋆ . The dependence of Z ⋆ (q i ) on the displacement along a particular vibrational mode, q i , is typically dominated by a linear term, which does not contribute to the renormalisation of Z ⋆ due to the symmetry of the vibrational wave function. Moreover, calculating Z ⋆ (q i , q j ) for simultaneous displacements q i , q j along pairs of vibrational modes suggests that cross-terms play a negligible role. Instead, the subdominant, higher-order terms in Z ⋆ (q i ) lead to a significant renormalisation of Z ⋆ , which is reflected in the increase in the nuclear polarisability by around 10% due to ZP nuclear vibrations. Evidence of this is supplied in the SI.
Vibrational corrections to the core electron density around the oxygen nuclei are too small to be resolved in the coarse-grained Fig. 4 due to the small root-meansquare displacements of the oxygen nuclei and affect neither the electronic quasi-particle band gap nor the polarisability of ice.
Vibrationally corrected permittivity. The experimentally observed permittivity, ǫ, of ice samples which contain various point and stacking defects, exhibits a nontrivial dispersion. The low-frequency behaviour around 50 Hz arises from proton-(dis-)order and Bjerrum and/or ionic defects. At intermediate frequencies between around 1 MHz and around 0.3 THz the relative permittivity is roughly constant at a value of around 3.2, and is well described as the sum of the the electronic response and the Debye relaxation of the dipoles of the water molecules. Above 0.3 THz crystal vibrational modes (or "pseudo-translations") begin to be excited, although the onset of the nuclear vibrational response occurs with the excitation of librational modes at above 10 THz. The nuclear vibrational response determines the detailed permittivity of ice up to around 100 THz before dying off, leaving only the electronic contribution to the permittivity of ice of around 1.7 arising from its electronic polarisability.
Here we focus on the permittivity of ice in the intermediate-frequency regime between 1 MHz and 10 THz, which we label simply as ǫ, although the electronic contribution to the permittivity provides insight into the optical permittivity of ice between 400 THz (red light) and 700 THz (blue light). The intermediatefrequency and optical regimes of the permittivity of ice play an important role, for example, in modelling cloud radiative processes and microwave propagation in atmospheric physics/climate science, as well as cloud measurements [45].
Recent work has probed both the permittivity of ice up to the THz regime [46] and its temperature dependence between 190 K and 260 K [47]. These experimental data provide reference points for computation and provide crucial information for understanding the role of proton-disorder in determining the permittivity of ice.
Here, the permittivities of the Monte Carlo sampled frozen-phonon configurations were calculated within the framework of Ref. [48] and using the Castep code to evaluate ⋆ and ǫ in the high-frequency limit, ǫ ∞ . Supplementary Fig. S5 shows ǫ for the 16 eightmolecule proton-orderings of Hirsch and Ojamae [39] and demonstrates that ǫ (unlike, for example, the electronic band gap [12]) depends sensitively on the particular proton-ordering. In contrast, the increase in ǫ due to nuclear motion is very similar for all proton-orderings and is larger than 5 % of the static-lattice value in all cases considered.
The temperature dependence of ǫ for a particular proton-ordering of Ih is negligible in comparison to the differences between values of ǫ of up to 18 % among proton-orderings.
The permittivity of a particular proton-ordering generally changes by no more than 0.5 % from zero temperature up to melting. This contrasts with the experimental increase in permittivity of protondisordered ice of close to 2 % over the same temperature range [46]. Boltzmann averaging over the 16 eightmolecule proton-orderings considered results in a temperature dependent ǫ which compares reasonably well with experiment (see Fig. 5), although it does not exhibit the clean quadratic increase of the experimental data. This is likely due to (a) the finite set of proton-orderings used [63] and (b) uncertainties in proton-ordering energetics. Nevertheless, Fig. 5 (a) highlights the importance of proton disorder and thermally induced changes in proton order.
In contrast to the total ǫ(= ǫ elec + ǫ dipole ), the electronic contribution to the relative permittivity, ǫ elec , of 1.87 (exp: 1.72 ± 0.02 [50]) is insensitive to the particular proton-ordering. ǫ elec differs by less then 0.3 % between different proton-orderings. The differences in the total ǫ arise from differences in the dipole contribution, ǫ dipole , between different proton-orderings. Also, ǫ elec has a negligibly small intrinsic temperature dependence, irrespective of the particular proton-ordering, so that ǫ elec exhibits a negligible temperature-dependence both for proton-ordered and proton-disordered ice.
Analogously to the mean ǫ, its anisotropy, ∆ǫ ≡ ǫ − ǫ ⊥ , where ǫ and ǫ ⊥ denote the permittivity along the principal axes parallel and orthogonal to the c-axis, varies substantially across proton-orderings and exhibits a negligible temperature dependence. Again, the Boltzmann average of the values of ∆ǫ compares reasonably well with the data for proton-disordered ice of Fujita et al. [47], as shown in Fig. 5 (b).
Both the static-lattice ǫ and its vibrational correction are sensitive to the choice of xc-functional. However, by constraining the simulation cell volume to the experimental value we largely eliminate the dependence on the xc-functional. The role of the xc-functional is further investigated in the Supplementary Information.

III. CONCLUSIONS
Proton-disorder, unlike temperature and pressure, substantially affects the proton position and momentum distributions in Ih. It is also crucial for the temperature dependence of the permittivity of Ih. Accurate staticlattice permittivities can be calculated for a particular proton-ordering using theoretical methods such as GW many-body perturbation theory [51], or quantum chemical methods, but the 5 % increase in the permittivity of Ih due to nuclear vibrational motion makes it imperative to account for vibrational effects when comparing to experiment.
The VSCF method provides a detailed understanding of the roles of anisotropy and vibrational anharmonicity for the proton momentum distribution, which will facilite the interpretation of experiments, for example, on nano-confined and supercooled water. Anharmonic vibrations are crucial in understanding the relative stability of Ih and Ic and many other phenomena. The VSCF vibrational method enables us to directly study their signatures in position and momentum distributions. Finally, the VSCF vibrational method provides the basis for efficient calculations of accurate position and momentum distributions and the means to understand the equilibrium proton dynamics in technologically important materials, such as KH 2 PO 4 [37] and other hydrogenbonded materials which are considered candidates for cheap, environmentally-friendly, organic electronics [52][53][54].
Acknowledgements. We acknowledge financial support from the Engineering and Physical Sciences Research Council (EPSRC) of the UK [EP/J017639/1]. The calculations were performed on the Cambridge High Performance Computing Service facility and the Archer facility of the UK's national high-performance computing service (for which access was obtained via the UKCP consortium [EP/P022596/1]).