Study of anharmonicity in Zirconium Hydrides using inelastic neutron scattering and ab-initio computer modeling

The anharmonic phenomena in Zirconium Hydrides and Deuterides, including {\epsilon}-ZrH2, {\gamma}-ZrH, and {\gamma}-ZrD, have been investigated from aspects of inelastic neutron scattering (INS) and lattice dynamics calculations within the framework of density functional theory (DFT). The observed multiple sharp peaks below harmonic multi-phonon bands in the experimental spectra of all three materials did not show up in the simulated INS spectra based on the harmonic approximation, indicating the existence of strong anharmonicity in those materials and the necessity of further explanations. We present a detailed study on the anharmonicity of zirconium hydrides/deuterides by exploring the 2D potential energy surface of hydrogen/deuterium atoms, and solving the corresponding 2D single-particle Schrodinger equation to get the eigenfrequencies. The obtained results well describe the experimental INS spectra and show harmonic behavior in the fundamental modes and strong anharmonicity at higher energies.

actinide hydride fuels [3]. Therefore, numerous experimental and theoretical investigations on their structures and phonon properties have been conducted and reported in publications since the 1950s.
The optical phonon properties of zirconium hydride were first analyzed using inelastic neutron scattering (INS) by Andresen et al. in 1957 [9]. Yamanaka [10] studied the thermal and mechanical properties of zirconium hydrides. Based on first-principles calculations, mechanical properties and thermodynamical properties of zirconium hydrides have also been broadly studied. Early studies by Ackland found that for ε-ZrH2, bistable crystal structures exist with minimal energy difference [11]. The subsequent studies found that one of the bistable structures may be unstable or metastable. Elsasser found that for the center of the Brillouin zone (Gamma point) of ZrH2, the potential energy profiles for hydrogen is a parabola up to about 300 meV and deviates at higher energies [12]. Besides theoretical studies, INS spectroscopy has also been employed to measure the vibrational spectra of ε-ZrH2. The split peaks of multi-phonon events were observed in 1983 by Ikeda [13], yet was not explained. Kolesnikov et al. later attributed those split peaks below multi-phonon events in γ-ZrH and γ-ZrD to the bound multi-phonon states [14,15].
Those experiments and calculations revealed that hydrogen vibrations in zirconium hydrides are highly anharmonic. Though phonon theory and techniques have been well developed, most of them are based on harmonic situations with a limited focus beyond that point. Phonon density of states (PDOS) with full consideration of anharmonicity can be extracted from molecular dynamics (MD) simulations by the velocity autocorrelation function (VACF); however, only frequencies information can be obtained in this case, which is insufficient for further phonon analysis. Recently, the temperature-dependent effective potential (TDEP) method has been developed to address the anharmonicity problem by the MD method and taking temperature effects into consideration [16,17]. However, these approaches are incapable of building directly comparable results with our broadband INS data for zirconium hydrides. Thus, an alternative technique is needed to address the problem. Besides, with a simple structure and large volume of experimental and simulation data, zirconium hydrides are excellent candidates to directly and unambiguously assess the quality of our anharmonic analysis.
Based on the considerations above, we present a qualitative and semiquantitative analysis that accounts for the effects observed in the INS spectra of zirconium hydrides by systematically investigating the potential energy surface (PES) of these materials using the density functional theory (DFT) and quantum theory. We sample PES in 2D planes due to considerations of the balance of practicality and accuracy. In this work, we successfully combined DFT and direct solutions of Schrodinger equations with INS spectra calculations. This work will hopefully improve our understanding of the mechanical and thermal properties of zirconium hydrides and can be applied to other metal hydride systems. The isotope effects were also investigated by checking the γ-ZrD system.

II. METHODS
Multiple ab initio simulation packages within the framework of DFT have been employed to carry out firstprinciples calculations in this work, including the Vienna ab initio simulation package (VASP) [18,19], and the Real-space Multi-Grid code (RMG) [20][21][22]. VASP is a plane-wave based DFT package with the implementation of multiple pseudopotentials. In all our VASP calculations, the electron-ion interaction and exchange-correlation functional were described by the projector-augmented wave (PAW) method [23] and the generalized gradient approximation (GGA) [24] with the parametrization of Perdew-Burke-Ernzerhof (PBE) [25], respectively. First-order Methfessel-Paxton scheme [26] with a smearing width of 0.05 eV was employed to integrate total energy in the Brillouin Zone (BZ), and the energy cutoff in the plane-wave functions was set to be 600 eV. Valence electron configuration in zirconium was carefully chosen as 4s 2 4p 6 4d 2 5s 2 to include semi-core states, and in hydrogen was 1s 1 . The self-consistent convergence of the total energy calculation is 10 -8 eV. RMG is an electronic structure calculation code based on the real-space method, implementing ultrasoft pseudopotential and norm-conserving pseudopotential. The RMG phonon methodology has been recently developed [27], and it is well suited for large-scale phonon calculations. In all RMG calculations, the real-space grid spacing was set to 0.15 Å, and norm-conserving pseudopotential with exchange-correlation described by PBE functional was used. Valence electrons for Zr were 4d 2 5s 2 . The convergence criteria for the root mean squared (RMS) change in the total potential energy per step was set as 10 -7 , and a Fermi Dirac occupation type with electron temperature as 0.05 eV was used in energy integration.

A. Atomic structures and phonon properties of ZrHx
The crystal unit cells of ZrHx and their phonon properties are shown in Fig. 1  respectively. The BZ paths and notation are adopted from [36]. corresponding to multi-phonon events excited from the ground state to the states higher than the first. The INS spectra simulated from VASP and RMG for ε-ZrH2 show excellent agreements with the experimental data at energies below 200 meV (see Fig. 3). However, discrepancies exist between the experimental data and the simulated INS spectra in the higher energy range. Overtone peaks at around 260 meV, 390 meV, and 510 meV exist only in the experimental data, suggesting the harmonic approximation is good for fundamental peaks but gets failing on higher energy excitations, and anharmonic effects need to be considered at higher energy. It should be noted that for γ-ZrH and γ-ZrD, similar split peaks were also observed in the experimental data (see Fig. 2).
The purpose of this work is to present a combined study regarding the anharmonicity of ZrHx by using firstprinciples calculations and neutron techniques.  we will try to simplify the complex many-body problem in two steps. By checking the structure of ε-ZrH2, we decided to sample the (112) plane (shifted so that the equilibrium position of the H atom is included in the plane) as this is one of the planes with the highest symmetry. The plane was sampled at a mesh density of 100 points per direction (mesh grid spacing around 0.005 Å, 10,000 mesh grids totally). Based on the sampled grid, a finer grid with a dimension of 1200×1600 = 1,920,000 was generated using a linear interpolation method to build the Hamiltonian matrix according to Eq. (2). A similar scheme has been applied to γ-ZrH and γ-ZrD. A view of the plane and the sampled PES is shown in Fig. 5. range of (0, 1) eV, and black lines are indicating the first six eigenfrequencies (see Fig. 6) from Schrodinger equations of ε-ZrH2 and γ-ZrD, respectively.
After solving the corresponding single-particle Schrodinger equation of H/D atom, the eigenfrequencies and wavefunctions which represent the energy states of atoms in the PES can be displayed, as plotted in Fig.   6. The three materials' eigenfrequencies have also been plotted in Fig. S2 [41] in the supplemental material, and direct comparisons between eigenfrequencies and experimental spectra can be found in Fig. 7. It can be seen that the anharmonic peaks in experimental spectra can be reasonably described by eigenfrequencies  Corresponding quantum numbers are denoted in [n, m], where n is in the Y direction and m is in the X direction. Notice that cross over happens between the quantum numbers of ε-ZrH2 and γ-ZrH, which are marked by orange arrows. The last two wavefunctions in γ-ZrD are degenerate. The Schrodinger equations were solved using [40]. reproduces the splitting of the hydrogen local mode peak at ~150 meV into two peaks (double and single degenerate, according to the symmetry of the H position), and the absence of dispersion due to disorder of the H atoms in the lattice (random defects). Notice that the relative higher intensity of low energy modes in experimental γ-ZrD0.99H0.01 than simulated low energy modes of γ-ZrD0.96875H0.03125 is mostly due to contributed low energy modes by α-Zr [14].
Note, that the exact knowledge of vibrational dynamics of different Zr-H materials is important not only from the point of view fundamental physics and chemistry but also for practical applications in the nuclear field such as critical safety studies. Therefore the obtained results on dynamics of ZrHx materials in the current study, which describes the anharmonicity of the INS spectra at energies above 300 meV, can be very useful for criticality and neutron transport simulations of the reactors containing Zr-H materials at anticipated operational temperatures up to 1200 °C (at which the first overtone level in Zr-H is ~10% populated).

IV. CONCLUSIONS
The anharmonicity phenomena in ε-ZrH2, γ-ZrH, and γ-ZrD have been identified and thoroughly studied with techniques from both DFT and INS. The anharmonicity is not apparent on the fundamental vibrations, and it becomes evident at higher energies. This effect is not appreciable around room temperature. However, ZrHx has been considered as a neutron moderator material for nuclear reactors; in this case, the anharmonic effects cannot be ignored since neutron with high energies will scatter from the hydrogen atoms and will experience deviations from the scattering kernel predictions based on the harmonic approximation. While V.