A Modeling and Neutron Diffraction Study of the High Temperature Properties of Sub-Stoichiometric Yttrium Hydride for Novel Moderator Applications

Low-enriched-uranium (LEU) reactor systems utilize moderators to improve neutron economy. Solid yttrium hydride is one of the primary moderator candidates for high-temperature (>700 °C) nuclear reactor applications. This is due to its ability to retain hydrogen at elevated temperatures compared to other metal hydrides. For reactor modeling purposes, both neutronic and thermos-mechanical modeling, several high-temperature properties for sub-stoichiometric yttrium hydride (YH2−x) are needed. In this paper, we present an atomistics and a neutron diffraction study of the high-temperature properties of Y  and  YH2−x. Specifically, we focus on the thermal lattice expansion effects in yttrium metal and yttrium hydride, which also govern bulk thermal expansion. Previously reported physical and mechanical properties for sub-stoichiometric yttrium hydride at ambient conditions are expanded using lattice dynamics to take into account high-temperature effects. Accordingly, an array of newly generated properties is presented that enables high-fidelity neutronics, and thermomechanical modeling. These properties include various elastic moduli, thermal expansion parameters for yttrium and yttrium hydride, and single-phase (YH2−x) and two-phase (Y + YH2−x) density as a function of stoichiometry and density.


Introduction
Microreactors are small, compact, and truck-transportable nuclear modular reactors. They are envisioned to produce between 1 kW and 10 MW of electrical power with selfregulating capabilities. This type of reactor is being considered for remote nuclear power production applications either as a hybrid source in conjunction with renewables or standalone. Due to their compact size and self-regulating characteristics, they are also considered for space-surface power and nuclear-electric propulsion applications. For civilian applications, low-enriched-uranium (LEU) fuel is required. LEU fuel necessitates the use of neutron moderators to reduce the size of the reactor core. For high-temperature (>700 • C) applications, the water moderators used for current nuclear power plants are not suitable and solid yttrium hydride is a main candidate as a primary moderator material due to its relatively low dissociation pressures at higher temperatures when compared to other metal hydrides. Some LANL microreactor designs that consider yttrium hydride as a moderating material include Opal, Zebra, Marvel, and Snowflake. A detailed discussion on the selection of this moderator material is provided by Mehta et al. [1]. Temperature gradients in the moderator result in the self-diffusion of hydrogen in metal hydrides due to Soret diffusion [2]. Accordingly, YH 2−x can transform into YH 2−x±∆ due to the hydrogen migration. As such, to model both the neutronic and mechanical behavior of the moderator in a reactor with high fidelity, when accounting for H diffusion effects, high-temperature stoichiometry-dependent properties for yttrium hydride are of paramount importance.
Research into using yttrium hydride as a nuclear reactor moderator started as a part of an aircraft nuclear propulsion program (ANP) [1]. In 1959, Lundin et al. [3] summarized a progress report on the yttrium hydride moderator material candidate for the ANP program. The thermal expansion data for yttrium hydride were included in this report from dilatometry, i.e., without considering possible phase transformations during heating or investigating the contributions from the hydrogen-charged metallic Y and YH 2 separately (as is needed for modeling). Under the same program, Parker [4] summarized the thermophysical properties of yttrium hydride. Following these efforts, in 1962 Flotow et al. [5] used adiabatic calorimeters to measure the heat capacities for YH 2 and YD 2 . In 1968, Mueller et al. [6] published the textbook Metal Hydrides. This book provided a comprehensive summary of fabrication, chemistry, testing, and material property measurements for all metal hydrides known at that time. Vetrano [7] reemphasized the use of hydrides as solid-state neutron moderators and reflectors for reactor applications. Stuhr et al. [2] studied the self-diffusion coefficient of yttrium hydride (1.80 ≤ H/Y ≤ 1.97) using neutron spectroscopy between temperatures 450 • C ≤ T ≤ 950 • C. H/Y represent the atomic ratio of hydrogen to yttrium in the yttrium hydride material. Yamanaka and group [8,9] used several experimental methods to study the mechanical and thermal properties of yttrium hydride close to stoichiometry (1.7 ≤ H/Y ≤ 2). Funston [10] and Beattie [11] used experimental methods to provide elastic properties of YH 1.5 and YH 1.93 , respectively. Funston's YH 1.5 Young's modulus data extended up to 1200 K. In the last few years, Yang et al. [12] and Schultz et al. [13] used atomistic simulation methods to study the mechanical properties of metal hydrides, including yttrium hydride. Peng et al. [14], Fu et al. [15], and Trofimov et al. [16] studied and summarized the thermodynamic information on the Y-H system, and thus the Y-H phase diagram. The work of Peng et al. [14] was used here to define the phase diagram for the Y-H system ( Figure 1). Recently, the combined experimental and theoretical work of Shivprasad et al. [17,18] extended the understanding of yttrium hydride fabrication, and the mechanical and thermophysical properties at H/Y ≥ 1.90. Mehta et al. [19] extended the stoichiometry-dependent structural and elastic properties at ambient conditions to 1.31 ≤ H/Y ≤ 2.0 using density functional theory (DFT) methods and neutron diffraction experiments. However, the thermal expansion effects were not included in that study, thus limiting that study to ambient conditions only. Hu et al. [20] also recently published a handbook of yttrium hydride material, which provides data for various moduli at ambient conditions. Mehta et al. [21] also studied the stoichiometry and temperature dependent single-phase heat capacity (c v ), and thermal scattering laws (TSLs) applicable for hightemperature nuclear reactor modeling. DFT quasi-harmonic approximation (QHA) for yttrium dihydride was first introduced here to validate the atomistic models using constant pressure heat capacity and thermal lattice strain.
In this work, the thermal expansion effects for stoichiometric yttrium hydride material are investigated and extrapolated for stoichiometries predicted to occur in the reactor application. Since yttrium hydride exists in both single-phase (at higher temperatures) and two-phase (at lower temperatures, see Figure 1) fields, for microreactor operation conditions, both the properties for single-phase yttrium metal (α-Y) and yttrium hydride (δ-YH 2−x ) are required to construct a comprehensive set of thermophysical properties. Accordingly, DFT-QHA simulations were performed for α-Y for the first time in this paper, and the δ-YH 2 DFT-QHA thermomechanical results are thoroughly analyzed, as they were excluded from our previous publication [20]. Furthermore, the results reported in our previous work [20] only contained the c p and the thermal strain data for δ-YH 2 but excluded various moduli as a function of temperature. Yttrium hydride samples were also investigated experimentally using neutron diffraction experiments and compared with the DFT results in both that paper and this one. Accordingly, thermophysical properties for yttrium metal (α-Y) are discussed in Section 3.1, yttrium hydride (δ-YH 2−x ) in Section 3.2, and two-phase material density in Section 3.3. All of the simulations presented here are part of an advanced multi-physics (DFT, neutronics, thermodynamics, and thermomechanics) model for nuclear reactor physics under development at Los Alamos National Laboratory.  [19,21] based on [14]). Crystal structures of α-Y and δ-YH 2 shown in respective phases. Large purple spheres are yttrium atoms, whereas small grey spheres are hydrogen atoms.

Methodology
The overarching atomic scale simulation method and experimental methods are described here. DFT methods were utilized for atomic scale modeling, whereas neutron diffraction experiments were utilized for validation of the DFT-generated properties. In this work, phonon density of states (DOS) for α-Y were evaluated using harmonic approximation, and the α-Y lattice expansion was evaluated using quasi-harmonic approximation. High-temperature neutron diffraction experiments were utilized to measure lattice expansion, site occupation factors, and atomic displacement parameters as a function of temperature of three different YH 2−x samples. Once the confidence in the DFT modeling methods were established for δ-YH 2 and α-Y by comparing to those experimental results, two-phase density for moderator applications was derived, which is crucial for reactor neutronics modeling and also serves as the most important result from this work.

Density Functional Theory
DFT is an efficient computational quantum mechanics method that predicts the electronic structure and interatomic forces of molecules and crystals, therefore enabling a comprehensive set of material properties for engineering applications. DFT is a widely used atomistics approach to determine the energy of many-body systems. In this study, we utilized a DFT-based quasi-harmonic approximation approach (QHA) to study yttrium metal (α-Y). We note that two additional volumes were added for yttrium dihydride (δ-YH 2 ) in this study with regard to our previous work [20], and the DFT-QHA analysis was reevaluated. The DFT code Vienna Ab initio Simulation Package (VASP) [22][23][24] was used to carry out the first-principles quantum mechanical simulations. The exchange correlational functional was described using the generalized gradient approximation (GGA) with Perdew-Burke-Ernzerhof (PBE) generated projected-augmented wave (PAW) potentials [25]. Partial occupancies for each orbital were set using first-order Methfessel-Paxton smearing [26]. For yttrium metal, the 4s, 4p, 5s, and 4d electron orbitals and the semi core s states were treated as valence. The collinear spin polarization was included for all calculations, and the projection operators were calculated in the reciprocal space to ensure the exact energy calculations within each electronic and ionic relaxation step.
An yttrium hydride (δ-YH 2 ) supercell of 2 × 2 × 2 fluorite unit cells, with 64 H and 32 Y atoms, was used in this study. A planewave cutoff energy of 500 eV with a 5 × 5 × 5 Monkhorst-pack k-point mesh within the Brillouin zone was utilized once the system energy was converged to 1 meV per atom. A Methfessel-Paxton smearing width of 0.2 eV was selected for all the QHA calculations. For energy minimization, the electronic loop and ionic loop were converged up to 1 × 10 −6 eV and 1 × 10 −5 eV, respectively. An yttrium metal (α-Y) hexagonal close-packed supercell of 3 × 3 × 3 unit cells with 54 atoms was chosen to model the two-phase region material density of sub-stoichiometric yttrium hydride. A planewave cutoff energy of 500 eV and 7 × 7 × 7 Monkhorst-pack k-point mesh with 0.2 eV Methfessel-Paxton smearing width was chosen after converging the system energy to 1 meV per atom. Similar to δ-YH 2 , the electronic loop and ionic loop were converged up to 1 × 10 −6 eV and 1 × 10 −5 eV, respectively. To predict the lattice expansion of stoichiometric YH 2 and Y, the quasi-harmonic approximation was utilized using the Phonopy code [27]. Using a stoichiometric δ-YH 2 supercell led to the creation of eleven supercells with distinct volumes using the VASP scaling parameter, adding two additional volumes beyond our previous work [20]. We observed that the addition of these two volumes had no impact on the coefficient of thermal expansion (CTE). For α-Y, seven distinct volume supercells were created. Using these different volumes, the harmonic approximation was applied to generate the phonon density of states through finite displacement perturbations for each volume. Once the phonon densities of states were produced for all volumes, the quasi-harmonic approximation was utilized to generate the thermal expansion parameters for both Y and YH 2 . The stoichiometry dependency of δ-YH 2−x at room temperature was adopted from our previous work [19] and reconstructed using DFT-QHA to account for the high-temperature effects, as discussed in this paper.

Neutron Diffraction
Yttrium hydride powders with compositions YH 1.48 , YH 1.81 , and YH 1.92 were prepared following procedures described in [17]. High-temperature neutron diffraction data was collected on the high-pressure/preferred orientation (HIPPO) neutron time-of-flight diffractometer [28,29] at the 800 MeV proton accelerator-driven short-pulsed neutron spallation source at the Los Alamos Neutron Science Center (LANSCE, Los Alamos, NM, USA) [30]. Details of the HIPPO instrument can be found in the previous references, and only a summary of the instrumental setup is provided here. The vanadium cans with the sample material were loaded in a vacuum furnace [31] with vanadium heating elements and heat shields, operating at a vacuum of~10 −4 Pa. One dataset per temperature, resulting from~3 h of count time (adjusted for proton current fluctuations from the linear proton accelerator), was analyzed. Neutron diffraction, similar to X-ray diffraction, allows for the extraction of information on, e.g., the phase fractions of the phases present in the sample (from overall scale factors of the peaks belonging to the individual phases), lattice parameters (from the peak positions), and unit cell content (from relative peak intensities).
The procedure applied to extract these parameters is as follows. The Rietveld method [32] as implemented in the General Structure Analysis System (GSAS) [33] was applied to analyze the data. Scripts developed in gsaslanguage [34] ensured that each sample at each temperature was analyzed with the same refinement approach. The initial structure for δ-YH 2 used for the refinement was the one reported for room temperature by Khatamian et al. [35] (space group F m −3 m, a = 5.197 Å), with both tetrahedral and octahedral sites defined. For α-Y, space group P 63/m m c and lattice parameters of a = 3.6474 A and c = 5.7306 A were used for the starting structure [36]. Thermal motion parameters U iso were constrained to be identical for both H sites but independent for all other atoms (Y atoms in YH 2 and α-Y) and initialized at a value of 0.01 A 2 . Diffraction data from the 140 • , 120 • , 90 • , and 60 • detector banks were refined simultaneously, excluding the lower resolution 40 • data. Leaving the YH 2 lattice parameter fixed at value of 5.2077 Å (deter-mined from ambient condition datasets with Cu as an internal standard, reported in [13]), the time-of-flight to d-spacing conversion factors were refined for all histograms, essentially recalibrating the sample position for each sample and providing reliable absolute lattice parameters at higher temperatures. Lattice parameters of YH 2 and α-Y, site occupation factors of tetrahedral and octahedral sites (without constraining the overall composition), and thermal motion parameters U iso for all atoms were refined. Absorption parameters and diffractometer constants were fixed at the values from the ambient condition refinement, and the ambient condition refinement was used as starting point for all other refinements with prescribed changes of the lattice parameters for temperatures above 700 • C. Examples of the refinements are shown in Figure 2. In Figure 2a, at room temperature (longer count time than the others) both α-Y and δ-YH 2 peaks are visible. In Figure 2b, collected at 800 • C, the dissolution of δ-YH 2 is observed as the α-Y peaks are stronger than δ-YH 2 peaks.

Results & Discussion
In the Section 3.1, thermal expansion in yttrium metal (α-Y) is described. Subsequently in Section 3.2, thermal expansion for δ-YH 2 is provided. Using the data from these two sections and previous work [20], in Section 3.3, temperature and stoichiometry-dependent single-/two-phase material density is provided.

Yttrium Metal
Displacement structures for α-Y were created for the VASP/PHONOPY toolkit. Under the harmonic approximation, the phonon density of states for yttrium were produced (shown in Figure 3) for each finite displaced structure. As expected, the phonons are in acoustic range for yttrium metal. The phonon density of states was then utilized to generate the thermodynamic properties such as constant-volume heat capacity and entropy ( Figure 4). Thermodynamic properties are in excellent agreement with the experiments by Jennings et al. [37]. The c v from this study is slightly overestimated (~3.5%) at~340 K when compared to the lattice heat capacity (c q ) from Jennings et al. [37], likely because c q ignores the electronic contribution. The entropies from this study and the literature were also compared to validate the DFT-generated α-Y properties, which are in excellent agreement with one another, also suggesting the validity of the generated phonons in this work.  After the DFT harmonic approximation was verified, the thermal expansion was investigated using the DFT-QHA method. DFT-QHA is necessary to generate lattice expansion-related properties, such as density and bulk modulus. Lattice expansions are of utmost importance for high-temperature applications. The physical density and CTE as a function of temperature were calculated and are provided in Figure 5. The density results were validated against the X-ray experiments reported by Spedding et al. [36]. Nevertheless, the density shown in the figure appears to have a large difference when compared to the experimental measurements, since the maximum difference is only~0.5%. All the validations provided here support the underlying DFT modeling utilized for α-Y metal for the purpose of this paper.

Yttrium Hydride
To understand temperature-dependent lattice expansion, earlier QHA models (similar to the ones utilized for α-Y) were adopted for YH 2 from our previous work [20]. In that work, only thermal lattice strain (ε) and constant-pressure heat capacity were provided strictly for validating the underlying phonons. However, in this work, remaining thermal expansion parameters including elastic moduli are extended and combined with those of α-Y in the two-phase field to enable the stoichiometry-and temperature-dependent understanding of the yttrium hydride material. Figure 6a presents the DFT-QHA-computed lattice parameters as a function of temperature for YH 2 . For validation, LANSCE neutron diffraction experiments are also shown. There is an over-prediction in the DFT lattice parameters of about 0.2% when comparing to the high-temperature neutron diffraction experiments. It is obvious from the neutron diffraction experiments that the lattice expan-sion of the YH 2 phase is stoichiometry-independent (stoichiometry of the entire compound, i.e., including varying amounts of α-Y). In other words, the CTE are independent of H/Y for yttrium hydride, as shown by Figure 6a. Since the CTE changes with varying stoichiometries are negligible, as shown by neutron diffraction experiments, and since QHA modeling requires extensive computational resources, QHA analysis was not extended for sub-stoichiometric yttrium hydride. As such, the CTE for YH 2 is utilized for all yttrium hydride sub-stoichiometries. The DFT-QHA approach also calculates the bulk moduli (B) as a function of temperature. The latter is essential for thermomechanical analysis for nuclear reactors. Since all overall stoichiometries share the same CTE for the δ-YH 2 phase up to H dissociation temperatures, the stoichiometry-dependent bulk modulus from Mehta et al. [19] in Figure 7a can be extrapolated to higher temperatures. In reality, the Poisson's ratio or ν will vary slightly as a function of temperature. However, if ν is treated as independent of temperature, Young's modulus or Y (Figure 7b) and shear modulus or G (Figure 7c) as a function of stoichiometry and temperature can also be obtained. This is possible due to the widely known relations described in Equations (1) and (2), below. The stoichiometry dependence of ν was also obtained from Mehta et al. [19].  We note that in the low-temperature two-phase region (Figure 1), the DFT-QHAgenerated single-phase moduli does not account for the two-phase grain interactions or two-phase CTE mismatch. This would require micromechanical modeling, which is outside the scope of this paper. Instead, it represents a hypothetical compound at low temperatures. However, the single-phase moduli are still applicable for high-temperature reactor applications when the material is in the single-phase field. All the moduli presented can be estimated using the equations provided below. The H/Y represents the overall ratio of hydrogen atoms to yttrium atoms and T is the temperature in K. These equations are valid for stoichiometries between 1.31 ≤ H/Y ≤ 2.0, and temperatures up to which hydrogen diffusion or dissociation (equilibrium pressure dependent) are observed. The room temperature errors with various publications in the literature are provided in Table 1. Note that the DFT-QHA data for various moduli are modeled in single-phase; however, the experimental data should exist in a mixture of binary phase according to the phase diagram in Figure 1. The DFT-QHA-model-predicted data is off by~6% and~15% on average compared to the most recent experimental data [17,20], as the experimental data contains a real two-phase system. The work of Schultz et al. [13] compares excellently with our data because of the size of their uncertainties. For the Young's modulus, our results are over-predicted compared to Funston et al. [10] and Hu et al. [20]; under-predicted compared to Setoyama et al. [8] and Beattie et al. [11]; and within the error bars for Shivprasad et al. [17] and Schultz et al. [13]. This shows that DFT-QHA results are within the acceptable range when compared to available experimental literature. The differences in our data when compared to the literature arise due to two primary assumptions made in our DFT-QHA models: (1) the single-phase single crystal ignores the mismatches from true binary phases since the α-Y is completely excluded, and its influence on binary-phase elastic properties are different than single-phase; (2) ν was treated as independent of temperature, which will slightly affect the DFT results between 0 K and 300 K. Typically, DFT and experimental results are compared at 0 K and room temperature, respectively; and (3) the utilization of the GGA-PBE method for DFT-QHA tends to under-predict the moduli. Even with discrepancies as high as 10% at room temperature, this data will prove to be useful at reactor temperatures (>~900 K) where the crystal is in single phase. Unfortunately, there is no experimental data to validate the high-temperature regime; however, this is where the DFT models from this study can be used to fill the gaps.

Two-Phase Density
The temperature-dependent density relations described in Figure 5 for α-Y and Figure 6 for δ-YH 2 , and the stoichiometry-dependent lattice parameter relation described by Mehta et al. [19], are used to derive the two-phase density (Figure 8). We emphasize that this result is different from Fig. 10 in [19], as in that paper the CTE effects were ignored. Yttrium hydride increases in density towards lower stoichiometry and lower temperatures in the two-phase region. At temperatures above 600 • C, hydrogen stoichiometries between 1.3 ≤ H/Y ≤ 2.0 are all in single phase. Single-phase YH 2−x increases in density with increasing H/Y and decreasing temperature. Combining the results from Mehta et al. [19] and Figure 6, the single-phase density for δ-YH 2−x or ρ δ−YH 2−x can be written in equation format, as shown in Equation (6). This equation is valid for stoichiometry between 1.31 ≤ H/Y ≤ 2.0 and temperatures up to which hydrogen diffusion or dissociation (equilibrium-pressure-dependent) are observed. The overall stoichiometry of the candidate moderator materials is around 1.8 with an expected gradient all within this validity range. The simplified equation below considers CTE effects due to temperature (in K) and stoichiometry effects due to the vacancy formation in single-phase yttrium hydride. The following equation also suggests that the density is more affected by temperature than the stoichiometry, ρ δ−YH 2−x g/cm 3 = 3.439 × 10 −11 · (T − 273.15) 3 −1.553 × 10 −7 · (T − 273.15) 2 − 1.555 × 10 −5 · (T −273.15) + 0.06047 · H/Y + 4.14096 (6) where T is in K and ρ δ−YH 2−x is in g/cm 3 . The two-phase density (ρ zY+(1−z)YH 2−x ) in equation format is more complicated as it contains the phase-diagram-dependent mixing parameter z, which is the weight fraction of yttrium metal (α − Y) in the two-phase system consisting of α − Y and δ − YH 2−x . The mixing parameter depends on the overall ratio of H/Y of the two-phase system, and the value of H/Y at the two-phase to single-phase boundary, i.e., (H/Y) PB . In addition, to derive density in the two-phase region, we assume that there are no H atoms in α − Y.
where T is in K and ρ zY+(1−z)YH 2−x is in g/cm 3 . z is given by The two equations provided can enable high-fidelity neutronics and criticality safety studies for the yttrium-hydride-moderated reactor, and is the most important result from this work. Without these equations, or the plotted density in Figure 8, a high-fidelity neutronics study is impractical. In addition, they provide an excellent basis for thermomechanical analysis.

Conclusions
For high-temperature reactor operations, yttrium hydride is being considered as a solid-state moderator. In order to successfully design and build such novel reactors, thermal expansion effects need to be accounted for. Accordingly, in this paper, temperaturedependent yttrium hydride properties are investigated and discussed. It was concluded that changes in hydrogen stoichiometry in yttrium hydride, as may occur in the reactor due to temperature gradients in the moderator, have no significant effect on the thermal expansion coefficients until high-temperature mechanisms related to the onset of the dissolution of the material are observed. Accordingly, the physical and structural properties of singlephase yttrium hydride are computed using DFT-QHA and compared with experiments in the literature. The density data were achieved at high confidence due to the mechanical properties data predicted here, since two-phase density was extrapolated using the QHA results for yttrium metal and yttrium hydride.