Benchmarking a ﬁrst-principles thermal neutron scattering law for water ice with a diffusion experiment

. The neutron scattering properties of water ice are of interest to the nuclear criticality safety community for the transport and storage of nuclear materials in cold environments. The common hexagonal phase ice I h has locally ordered, but globally disordered, H 2 O molecular orientations. A 96-molecule supercell is modeled using the VASP ab initio density functional theory code and PHONON lattice dynamics code to calculate the phonon vibrational spectra of H and O in ice I h . These spectra are supplied to the LEAPR module of the NJOY2012 nuclear data processing code to generate thermal neutron scattering laws for H and O in ice I h in the incoherent approximation. The predicted vibrational spectra are optimized to be representative of the globally averaged ice I h structure by comparing theoretically calculated and experimentally measured total cross sections and inelastic neutron scattering spectra. The resulting scattering kernel is then supplied to the MC21 Monte Carlo transport code to calculate time eigenvalues for the fundamental mode decay in ice cylinders at various temperatures. Results are compared to experimental ﬂux decay measurements for a pulsed-neutron die-away diffusion benchmark.


Introduction
Water is the most common neutron moderator in nuclear reactor applications. The dynamics and thermal neutron scattering properties of liquid water have been extensively studied for decades. Criticality safety analyses routinely examine the reactivity effect of a liquid water environment during the transport and storage of fissile material. For cold weather conditions, it is of interest to have a thermal neutron scattering law (TSL) for water ice which characterizes the impact of a chemically bound crystalline structure on neutron scattering and reflection.
In H 2 O molecules, four unbonded oxygen valence electrons form a region of negative charge opposite the two strong internal O-H covalent bonds. In condensed matter, H 2 O molecules participate in hydrogen bonding with each other due to electrostatic dipole-dipole interactions [1]. These relatively weak hydrogen bonds are responsible for small variations in the H-O-H bond angle for different phases and temperatures. In the solid form, hydrogen bonds are also responsible for the geometric flexibility of H 2 O alignments in various phases of ice [2]. The details of the particular interatomic structure affect the available vibrational modes with which thermal neutrons may interact.

The structure and modeling of ice Ih
Hexagonal ice Ih is the phase of water ice that exists under common temperature and pressure conditions. While the oxygens lie on a hexagonal lattice, the ice Ih structure exhibits locally ordered but globally disordered H 2 O molecular alignments [3]. Thus, ice Ih may be considered a e-mail: jesse.holmes@unnpp.gov to have an "average" structure which differs from the local structure [4].
There have been long-standing challenges associated with determining the vibrational phonon energy spectrum for ice Ih due to its globally varying structure. Since H 2 O has a very large incoherent to coherent nuclear scattering cross section (XS) ratio, inelastic neutron scattering can be an effective means of experimentally probing spectral features associated with the average structure, but the full one-phonon spectra for H and O (which are required for TSL calculations) can be difficult to extract accurately [5,6].
Lattice dynamics models have been used for decades to predict the vibrational modes of materials based on force constants. Historically, including for ice Ih, these force constants were fitted to thermodynamic measurements and other experimental data [5,7,8].
In recent years, increased computational power has made the use of first-principles density functional theory (DFT) attractive for investigating the electronic and structural properties of crystalline solids. While such studies for ice Ih have been complicated by the difficulty in modeling the hydrogen bond dipole-dipole interactions and the disordered hydrogen geometry, several published results have compared well to experimental inelastic scattering spectra [9][10][11]. Within this methodology, the phonon spectra can be determined by solving a dynamical matrix for calculated interatomic force constants [12,13].
In the present work, the phonon spectra for H and O in ice Ih are calculated for the primary purpose of producing an accurate thermal neutron scattering kernel [14] for ice Ih. In 2002, Bernnat et al. [15] produced a TSL for this structure using a smoothed representation of the phonon spectrum generated by Nakahara [8] in 1968 by a historical  root sampling method. However, there is no published Evaluated Nuclear Data File (ENDF) [16] from this work, and there are no published works generating a high-quality thermal neutron scattering kernel for ice Ih from modern first-principles methods.
The structure of ice Ih is modeled using the VASP ab initio DFT code [17][18][19] with a 12-molecule unit cell as shown in Fig. 1. A 5 × 5 × 5 k-mesh and electron plane-wave energy cutoff of 900 eV is used to relax the system. The generalized gradient approximation (GGA) is used for the electron exchange-correlation energy, which is known to be significantly superior to the local density approximation (LDA) in predicting hydrogen bond lengths. Even so, GGA and the available atomic potentials for H and O significantly underestimate hydrogen bond lengths. This is a known deficiency of DFT [9]. To compensate for this, high-accuracy lattice constants for ice Ih from Röttger et al. [20] at 115 K are applied, constraining the unit cell dimensions. Since the thermal expansion of ice Ih is negligible from 0 K to 115 K (the temperature at which the Torres et al. [21] experimental total cross sections are available), this procedure is justified.
Interatomic forces are calculated using a 96-molecule 2 × 2 × 2 supercell on a 1 × 1 × 1 k-mesh with a 400 eV plane-wave cutoff energy. Using supplied force constants, the lattice dynamics code PHONON [22] solves eigenvalues of the dynamical matrix for randomly sampled k-points in the first Brillouin zone, collects the results, Figure 3. Experimentally measured inelastic neutron scattering spectrum for ice Ih [6]. The x-axis energy scale differs from Fig. 2. The y-axis units are arbitrary. and generates the required phonon spectra. In Fig. 2, the phonon spectra of H and O are shown together, where H is normalized to 2/3 and O is normalized to 1/3. Note the distinct translational, librational, bending and stretching modes as energy increases. Figure 3 compares the calculated results to experimentally measured inelastic neutron scattering data for low energy transfers [6]. While the calculated features are correctly positioned in energy (the energy scales differ in the two figures), the librational region is clearly broadened for the experimentally measured average structure [9]. The calculated high-energy bending and stretching modes match experimentally measured peak locations and are largely decoupled from molecular alignments. [23,24]. Figure 4 shows three theoretical total XS calculation comparisons to the experimentally measured total XS for ice Ih at 115 K [21]. The calculated phonon spectra were used to determine inelastic and elastic scattering in the incoherent approximation (since H 2 O is largely an incoherent scatterer) using NJOY2012 [14]. Next, absorption was added in [25] to give the total XS. The C/E results are shown in Fig. 5. Optimizing the weight of the librational region to account for an average structure improved consistency with total experimental XS data from the light blue (top) line to the red (middle) line.

Comparison of theoretical and experimental total cross sections
Elastic scattering by oxygens is physically entirely coherent. The incoherent approximation overestimates oxygen elastic scattering at lower energies. However, this model is considered acceptable since the oxygen scattering XS is very small relative to hydrogen. While the removal of oxygen elastic scattering (in the incoherent approximation) appears to improve results further -dark green (lower) total XS curve in Fig. 4 -the true total cross sections may fall somewhere between this and the red (middle) curve. As the red curve scattering kernel yielded excellent results in the diffusion benchmark tests (better than the dark green curve scattering kernel), it is considered the final optimized scattering kernel for ice Ih. In any case, considering that the quoted ENDF data uncertainty on the free scattering cross section of H 2 O exceeds 2% [25], and the experimental data in Fig. 4 is normalized to this quantity, the XS results are reasonable.

Comparison of theoretical and experimental time eigenvalues
A pulsed neutron die-away experiment was used to study the temperature dependence of neutron diffusion in ice Ih  by measuring the fundamental spatial mode decay time eigenvalues for cylinders of ice of various sizes [26,27]. The neutron population persisting in a system after an initial pulse will equilibrate into the fundamental spatial mode as the population decays. In the one-speed diffusion model, the flux will be of the form where a is the macroscopic absorption cross section D = 1/ transport is the diffusion coefficient, v = the effective average neutron velocity, B 2 = geometric buckling, and C = the diffusion cooling coefficient.
Geometric buckling is a quantification of the curvature of the fundamental mode shape. High bucklings are associated with small geometries and greater leakage rates. For H 2 O, thermal scattering cross sections are much larger than absorption cross sections. Therefore, the v a pure absorption term will only be dominant for verylow buckling cases (associated with large geometries),  From top down, the green line is a free gas library at ice density, the brown line is a subcooled liquid water library [25] at ice density, the red line is the final optimized library, the black line is experimental data [26], and the blue line is the originally generated theoretical library without optimization.
where the thermal scattering kernel applied will be relatively unimportant. When buckling is high (for small geometries), the thermal scattering kernel has a dominant role in determining the time eigenvalue, defined as α = v a + vDB 2 − CB 4 . The diffusion cooling term captures the geometric effect that higher energy neutrons will be preferentially lost at the leakage boundaries, resulting in an average neutron population at a temperature lower than thermodynamic equilibrium. ENDF File 7 [14,16] thermal scattering libraries for ice Ih are generated using the THERMR module of NJOY2012 [14]. These are processed by NDEX [28] and supplied to the MC21 transport code [28], in which the ice cylinders are modeled. Calculated time eigenvalues for ice Ih, as well as for free gas and subcooled liquid water [25] models, are compared to experiment [26]. See Table 1 and Figs. 6 and 7 for detailed results data.  [26]. Units of B2 are in cm −2 . Units of α are in sec −1 . Uncertainties are in absolute units.

Conclusions
A thermal neutron scattering law for ice Ih has been generated through ab initio methods and optimized incorporating experimental lattice constants, experimental inelastic neutron scattering spectra, and experimental total cross sections. The final resulting scattering kernel performs very well in reproducing experimental measurements of fundamental mode decay time eigenvalues for a pulsed-neutron die-away diffusion benchmark for a variety of temperatures and geometries. The average C/E over all temperatures for the six highest buckling geometries was 1.0032. Liquid water and free gas libraries perform much more poorly, as seen in Fig. 7. The reflective and scattering properties of the complicated ice Ih structure with globally disordered molecular alignments have been sufficiently well characterized to predict a wide variety of physical phenomena which depend on the interaction of thermal neutrons with the structure's vibrational modes.