Generation of the TSL for Zirconium Hydrides from Ab Initio Methods

Zirconium hydride (ZrHx) is a moderator material used in TRIGA and other reactors that may exist in multiple phases with varying stoichiometry, which include the δ phase and the ϵ phase. Current ENDF/B-VIII.0 ZrHx thermal scattering law (TSL) evaluations do not distinguish between phases. These sub-libraries were generated with the LEAPR module of NJOY using historic phonon spectra derived from a central force model and assume incoherent elastic scattering for both bound hydrogen and zirconium, which neglects the effects of crystal structures important for scattering from zirconium bound in ZrHx. In this work, the TSLs for hydrogen and zirconium bound in δ-ZrHx and ϵ-ZrH2 were generated from phonon spectra derived from modern ab initio lattice dynamics methods and ab initio molecular dynamics. Subsequently, TSLs for hydrogen and zirconium in ZrHx and ZrH2 were generated using the Full Law Analysis Scattering System Hub (FLASSH) code. The built-in generalized coherent elastic routine was used to generate the previously neglected elastic contribution from zirconium for this material. The present TSLs provide both a re-evaluation of the current ZrH sub-libraries and expansion of the set of TSLs available for the examination of neutrons in systems with zirconium hydride, permitting explicit treatment of δ and ϵ phases.


Introduction
Zirconium hydride is used as a moderator material in TRIGA reactors and has been historically used as a moderator in other reactors requiring high proton density for moderation (e.g., SNAP-10A [1]). Additionally, ZrH x is a corrosion product in Zircaloy. ZrH x exists in multiple phases or polymorphs with varying stoichiometry [2]. The most significant of these phases are the δ phase and phase, the former of which is the primary constituent of TRIGA fuel and Zircaloy corrosion [3,4]. At room temperature, the δ phase is dominant for 1.59 < x < 1.64 and exists in a mixture with α-Zirconium for x < 1.59 below the β-Zr ↔ α-Zr + δ-ZrH x reaction eutectoid temperature of approximately 820 K [2]. The phase is dominant for x > 1.74 and is mixed with the δ phase for 1.64 < x < 1.74. This stability range for the phase is approximately from room temperature up to nearly 800 K; however, the stability range for the δ phase grows considerably with increasing temperature, shrinking the mixed-phase region [2]. As a result of crystal binding of hydrogen in ZrH x , energy exchange associated with thermal neutron scattering involves narrow, well-separated energy levels in contrast to the large continuum spectra observed in other materials [5,6]. This behavior is characteristic of metal-hydrides [5,6] and greatly influences neutron thermalization effects, such as the warm neutron principle responsible for the large negative moderator feedback coefficient of TRIGA reactors.
In neutron transport calculations (e.g., Monte Carlo) of reactors, Evaluated Nuclear Data File (ENDF) thermal scattering law (TSL) sub-libraries are used to capture the effects of crystal binding effects on low-energy neutron scattering. The TSL, S(α, β), is a temperaturedependent probability distribution defining the possible scattered states for an incident neutron, and has contributions from distinct effects (d) due to scattered wave interference between different nuclear scattering sites and self (s) non-interference effects [7], The momentum and energy transfer resulting from scattering are described by the unitless parameters α and β as, with µ and A representing the scattering cosine and nuclide to neutron mass ratio, respectively. The incident and scattered neutron energies are E and E , and k B and T are the Boltzmann constant and temperature. In crystals, the phonon expansion is used to compute the TSL such that [7,8], which is a summation of increasing orders of the convolution of the vibrational (phonon) spectra, with p = 0 representing elastic scattering and p > 0 representing inelastic scattering. In current ENDF thermal scattering sub-libraries, the incoherent approximation is applied, such that distinct effects on inelastic scattering are neglected. Within this approximation, the double-differential scattering cross section is related to the TSL using Fermi's Golden rule as [7], where σ coh and σ inc are the bound coherent and incoherent nuclear cross sections of the nuclide, and T is the temperature of the material. The current ENDF/B-VIII.0 ZrH x TSLs were generated in the incoherent approximation within the LEAPR module of NJOY using the vibrational spectra calculated with the central force model developed at General Atomics (GA) [8][9][10], which is based upon the phase (i.e., ZrH 2 ). Nevertheless, these TSLs have historically been used to represent both phases in reactor physics calculations. Elastic scattering was assumed to be incoherent for both hydrogen and zirconium, such that S 0 d = 0; an approximation which is reasonable for hydrogen due to its large incoherent nuclear scattering cross section. This approximation, however, neglects the effects of crystal structure important for scattering from zirconium bound in ZrH x .
In this work, the TSLs for hydrogen and zirconium bound in ZrH x are generated from phonon spectra derived from ab initio lattice dynamics (AILD) and ab initio molecular dynamics (AIMD) for both δ-ZrH x and -ZrH 2 , thereby expanding the set of available sub-libraries for this material. These TSLs were generated using the Full Law Analysis Scattering System Hub (FLASSH), which includes a generalized coherent elastic scattering capability [11]. At this stage, coherent elastic scattering only incorporates zirconium, enabling the sub-libraries to be extensible, within the limitations of the ENDF format, to the entire stoichiometry domain of each phase and mixed phase.

Computational Methods
The δ and phases of ZrH x are characterized by cubic and tetragonal crystal structures, respectively [12]. For the purposes of ab initio calculations, the stoichiometry is represented as an idealized x = 1.5 for δ phase and x = 2 for phase. The unit cells of each phase used in simulation are shown in Figure 1. Initially, the geometry of these structures were optimized through electronic structure calculations using the Vienna ab initio Simulation Package (VASP) [13,14] with plane-augmented wave (PAW) potentials and a GGA-PBE exchange correlation functional [15,16]. A plane wave cutoff of 450 eV and k-mesh of less than 1.8 nm −1 were used and found to be sufficient for energy convergence of less than 4 meV/atom. The predicted lattice parameter of the δ phase are a = 0.477 nm, which is in good agreement with the experimental value of a = 0.478 nm [12]. The predicted phase parameters of a = 0.354 nm and c = 0.440 nm are also in reasonable agreement with experiment, a = 0.352 nm and c = 0.445 nm [12]. The phonon density of states (DOS), which is the fundamental input to the generation of TSL in the incoherent approximation, was generated in the harmonic approximation for δ-ZrH 1.5 and -ZrH 2 using AILD methods involving the coupling of the VASP [13,14] and PHONON [17] codes, described in Ref. [18]. These calculations were performed with the MedeA simulation environment [19]. Hellmann-Feynman forces were extracted from VASP calculations of 2 × 2 × 2 (80 atoms) δ-ZrH 1.5 and 3 × 3 × 2 (108 atoms) -ZrH 2 supercells with the previously listed convergence parameters. The forces describing the small displacements of zirconium atoms are expected to be well-characterized within the harmonic approximation; however, hydrogen-due to its small mass and large vibrational displacement-may have significant anharmonic contributions to its motion even at 0 K (i.e., zero-point motion). AIMD simulations, which inherently include anharmonicity, were therefore performed using methods developed in previous work [20,21] to generate the hydrogen DOS. Using a reduced 350 eV planewave cutoff and a 1 × 1 × 1 k-mesh, a 3 × 3 × 3 (270 atoms) δ-ZrH 1.5 and 4 × 4 × 3 (288 atoms) -ZrH 2 supercell were first allowed to thermodynamically equilibrate for 1 ps at 0.00025 ps time-steps under a 300 K NVT thermostat. Atomic trajectories were then evolved for an additional 1.5 ps. Subsequently, the velocities at each time-step during this second time-period were used to generate the phonon DOS for hydrogen as the Fourier transform of velocity auto-correlation function, which, at 300 K, captures intrinsic anharmonic effects from the atomic forcefield without introducing significant additional spectral broadening.
The resulting hydrogen phonon DOS from AIMD, as seen in Figure 2, are shifted to lower energies compared to AILD calculations, and are more consistent with inelastic neutron spectroscopy (INS) measurements. Additionally, the AIMD DOS are more analogous to experiments than the GA spectrum used in the current ENDF evaluation [8,10], which has been previously noted to be qualitatively different from the experimental spectra [5,22]. Although the difference between the approximate AILD and AIMD average phonon energies is less than 0.01 eV, the deviation has an increasingly large impact on the predicted TSL and corresponding cross section for neutron scattering events with hydrogen involving the emission of multiple phonons. Specifically, the predicted oscillations using an AILD spectra diverge from experimental observation with increasing phonon order. Consequently, the AIMD phonons for hydrogen are expected to be more suitable for use in the generation of TSL.  [5]. The experimental, AILD, and AIMD phonon spectra differ between phases. The AIMD spectra is used for subsequent evaluation of the hydrogen TSLs in this work.
TSL for zirconium and hydrogen in δ-ZrH x (Zr(ZrH x ) and H(ZrH x )), as well as zirconium and hydrogen in -ZrH 2 (Zr(ZrH 2 ) and H(ZrH 2 )) were generated in the incoherent approximation as File 7 (MF = 7) thermal scattering sub-libraries (NSUB = 12) in ENDF-6 format [23] using the FLASSH code [11]. Each File 7 included evaluation of the TSL for temperatures between 296 K and 1200 K. The bound and free nuclear scattering cross sections (1 b = 10 −28 m 2 ) used for each element are listed in Table 1 and were derived from relevant data in the NIST neutron scattering length database and ENDF/B-VIII.0 neutron reaction sub-libraries (at E = 0.0253 eV from File 3 MT = 2 processed at 0 K) [24,25]. The zirconium File 7 corresponds to the natural isotopic mixture, as reported in IUPAC atomic weight and isotopic abundance evaluations [26,27]. The phonon DOS used for zirconium and hydrogen were from AILD and AIMD simulations, respectively. A phonon order of 300 was used in the phonon expansion of all inelastic reactions (MT = 4). Table 1. Elemental scattering cross sections used in the present TSL evaluations for ZrH x . The coherent and incoherent cross sections are bound. Each free atom cross section represents an isotopic average of the ENDF/B-VIII.0 neutron reaction sub-libraries at 0.0253 eV [25]. For hydrogen, the incoherent cross section corresponds to the Sears database [24], whereas the coherent scattering cross section is estimated as the difference between the bound total and bound incoherent scattering cross sections. The coherent elastic scattering due to crystalline structure of each phase was computed using the generalized elastic scattering treatment in FLASSH and stored on the corresponding elastic reaction (MT = 2) of the zirconium File 7. The Debye-Waller matrix in the present work was computed in the cubic approximation. Although hydrogen has an appreciable coherent scattering length (−3.7 fm) when compared to zirconium (7.3 fm) [24,25], the contribution of the hydrogen sublattice was neglected so that the TSL is extensible across the entire stoichiometric range of each phase.
The hydrogen TSLs used the incoherent nuclear scattering option for MT = 2. The nuclear scattering cross section stored on this reaction is only the nuclear incoherent scattering cross section, whereas the total free-atom cross section is used for MT = 4 per Equation (5). This treatment, while in contrast to the use of the total cross section for MT = 2 in ENDF/B-VIII.0, is both consistent with the neglection of the hydrogen contribution to coherent elastic scattering and the treatment of nitrogen in uranium mononitride [25], which like hydrogen, has a non-negligible coherent and incoherent nuclear scattering cross section.

Thermal Scattering Laws
Inelastic neutron scattering techniques provide a direct measure of the double-differential cross section of a compound, which is also directly calculable from the TSL via Equation (5). The double-differential scattering cross section for ZrH 2 generated in the present work at the Naval Nuclear Laboratory is compared to the experiment and ENDF/B-VIII.0 in Figure 3. The anticipated oscillatory behavior due to the localized, higher-energy hydrogen phonon DOS is observed, whereby peaks occur at intervals of approximately 0.14 eV of energy loss. Both the Naval Nuclear Laboratory and current ENDF TSLs demonstrate reasonable agreement with experiment [28]. Despite the similarity in the energy structure of the present and ENDF/B-VIII.0 double-differential cross sections, a slight divergence occurs in each subsequently increasing energy loss peak due to differences in the average hydrogen phonon energy; however, as noted previously, the Naval Nuclear Laboratory phonon DOS more closely resembles the experiment.
The differential scattering cross section of -ZrH 1.92 , shown in Figure 3, is computed as an integral of Equation (5) over all scattered neutron energies. The Naval Nuclear Laboratory TSL and ENDF TSL provide similar predictions and trend well with experiment [29]. These calculations exclude the contribution of coherent elastic scattering from the differential cross section; therefore, the locations and relative amplitudes of the coherent elastic reaction (i.e., diffraction pattern) of the Naval Nuclear Laboratory zirconium TSL are also shown. Deviation from the experiment in the calculated differential cross sections at low scattering cosine are due to the use of the incoherent nuclear scattering cross section as opposed to the total elastic reaction, MT = 2, of the Naval Nuclear Laboratory TSLs. At low-incident neutron energies where coherent elastic scattering is significant, diffraction peaks occur in the experimental data that result in a deviation from the continuous differential behavior of the calculations. Overall, the calculated diffraction pattern matches these peaks, thereby demonstrating the importance of this contribution for the prediction of angular distributions.
The inelastic contributions to the Naval Nuclear Laboratory TSLs for both δ-ZrH x (H(ZrH x ), Zr(ZrH x )) and -ZrH 2 (H(ZrH 2 ), and Zr(ZrH 2 )) are illustrated in Figure 3 at temperatures of 296 K and 700 K, where the latter is representative of the maximum fuel temperature of a 2000 MW pulse in a TRIGA reactor [3]. In general, the TSLs for hydrogen and zirconium are similar between the two phases. The zirconium TSLs indicate a quasi-elastic scattering spectrum at low momentum transfers, trending toward free gas behavior at large momentum transfer (e.g., α = 1.979). Oscillations in the hydrogen TSLs strongly favor energy loss (β < 0) at 296 K, but have a significant contribution from 0.14 eV 1-phonon up-scattering events at 700 K. This behavior is a contributor to the strong negative moderator feedback coefficient of TRIGA fuel systems. At low momentum transfers (e.g., α = 0.4943), 1-phonon events are dominant; however, as the momentum transfer increases the possibility of down-scattering, involving multiple phonons becomes likely, exemplified by α = 5.067. However, the multi-phonon peaks occurring at integer values of approximately 0.14 eV are similar for the two phases, and oscillations diverge with increasing |β| due to the small increase in the average phonon energy of δ relative to the phase, that is also observed in the experiment [5]. . Computed double-differential cross section of -ZrH 2 for E = 0.572 eV and µ=0.766 (top left) and differential cross section of -ZrH 1.92 for E = 0.0225 eV (top right) compared to the experiment, Purohit [28] and Korbichler [29]. TSL for H (bottom left), and Zr (bottom right) in δ-ZrH x and -ZrH 2 are shown at 296 K and 700 K. Both α and β are referenced to 293.6 K (i.e., k B T = 0.0253 eV). The double-differential scattering cross sections for the present Naval Nuclear Laboratory (NNL) and ENDF/B-VIII.0 TSL evaluations were generated with the FLASSH code [11]. The diffraction peaks in the differential scattering cross section indicate the location of contributions from coherent elastic scattering in the NNL evaluation that are not included in the NNL plot. Diffraction presents as a deviation from continuous behavior in the experimental differential cross section. Experimental double-differential cross sections were obtained from the time of flight measurements. Integrated total scattering cross sections for δ-ZrH x and -ZrH 2 were generated using the NDEX nuclear data processing package [30] and are shown in Figure 4. NDEX uses the β grid tabulated on the TSL and an adaptive energy mesh to provide high-resolution cross sections [31], which is important for predicting the oscillator behavior of metal hydride moderators that is not fully captured with the coarse, fixed energy meshes used in other codes (e.g., NJOY [9]), as demonstrated in previous work [31,32]. Several compositions were selected for comparison between the two models and experiment, with x = 1.6 for the δ phase and x = 2 for phase representing TRIGA and SNAP-10A fuel, respectively [1,3]. Each cross section is illustrated at room and elevated temperatures, with 700 K corresponding to pulse characteristics of a TRIGA and 800 K corresponding to the approximate maximum operating temperature of SNAP-10A [1].
At 296 K, the impact of changes in the coherent and incoherent elastic scattering is appreciable for both phases, resulting in a 6-10 b difference in the total scattering cross section below the Bragg edge between the Naval Nuclear Laboratory and current ENDF evaluations. In particular, the Bragg peaks at 0.0027 and 0.0036 eV in the δ-ZrH Naval Nuclear Laboratory evaluation are in reasonable agreement with experimental data [3], but are absent from the current ENDF evaluation. Below the Bragg cut-off (0.0027 eV), the experimental cross section is in between the Naval Nuclear Laboratory and ENDF calculations, consistent with the neglect of coherent elastic effects of hydrogen (e.g., correlated hydrogen vacant sites). Additionally, ZrH 1.6 shows a parallel 1/v trend between the Naval Nuclear Laboratory and ENDF evaluations in this region, whereas ZrH 2 diverges with less pronounced 1/v behavior in the Naval Nuclear Laboratory evaluation. This disparity arises from a decreased contribution of low-energy (<0.03 eV) phonons in the phase compared to the δ phase, present in both AILD and AIMD calculations.
As temperature increases, inelastic scattering with hydrogen begins to dominate, and the 0.14 eV 1-phonon up-scatter becomes the primary contributor to the 1/v scattering cross section. Consequently, the total scattering cross sections for the Naval Nuclear Laboratory and ENDF evaluations are nearly the same at both 700 K and 800 K; however, an observable consequence of coherent elastic scattering on the total scattering cross section persists between approximately 0.001-0.004 eV. At incident energies above 0.05 eV, both calculated cross sections trend well with the experiment. In this range, inelastic scattering causes oscillations in the cross section at intervals of approximately 0.14 eV that broaden with increasing temperature, and quantitative differences between Naval Nuclear Laboratory and ENDF evaluations arise primarily from a 0.4 b difference in the free-atom cross section used in the generation of the zirconium TSLs. . Scattering cross section of δ-ZrH compared to experiments at room temperature, Whittemore [6] and Schmidt [33] (left). Scattering cross sections of δ-ZrH 1.6 and -ZrH 2 compared to current ENDF/B-VIII.0 evaluations (right). The integrated cross section for both the NNL TSL evaluation and ENDF/B-VIII.0 were generated with the NDEX nuclear data processing code [30,31].

Conclusions
Phonon spectra for hydrogen and zirconium in δ and zirconium hydride were developed using predictive AIMD and AILD simulations. Subsequently, these spectra were used in the FLASSH code to generate new TSL evaluations for zirconium hydride which include, for the first time, the previously neglected contribution of coherent elastic scattering in multiple phases of this material. These evaluations include both δ-ZrH x (H(ZrH x ) and Zr(ZrH x )) and -ZrH 2 (H(ZrH 2 ) and Zr(ZrH 2 )), permitting explicit modeling of each phase in reactor physics calculations or benchmarks (e.g., TRIGA, SNAP-10A) across the entire range of their stoichiometry. Currently, only a single sub-library set (H(ZrH) and Zr(ZrH)) is included in the ENDF evaluation for this material, such that the new Naval Nuclear Laboratory evaluations may be considered not only a re-evaluation, but also an expansion of the ENDF/B library database. Moreover, differences in both the coherent elastic reaction and hydrogen phonon spectra between the evaluations provides the capability to study the importance of separate treatment for the δ and phases in neutronics calculations.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.