Study of Anharmonicity in Zirconium Hydrides Using Inelastic Neutron Scattering and Ab-Initio Computer Modeling

Abstract

The anharmonic phonon behavior in zirconium hydrides and deuterides, including ϵ-ZrH₂, γ-ZrH, and γ-ZrD, has been investigated from aspects of inelastic neutron scattering (INS) and lattice dynamics calculations within the framework of density functional theory (DFT). The harmonic model failed to reproduce the spectral features observed in the experimental data, indicating the existence of anharmonicity in those materials and the necessity of further explanations. Here, we present a detailed study on the anharmonicity in zirconium hydrides/deuterides by exploring the 2D potential energy surface of hydrogen/deuterium atoms and solving the corresponding 2D single-particle Schrödinger equation to obtain the eigenfrequencies, which are then convoluted with the instrument resolution. The convoluted INS spectra qualitatively describe the anharmonic peaks in the experimental INS spectra and demonstrate that the anharmonicity originates from the deviations of hydrogen potentials from quadratic behavior in certain directions; the effects are apparent for the higher-order excited vibrational states, but small for the ground and first excited states.

1. Introduction

Metal hydrides have been studied extensively during the past decades. In particular, there are special interests in their properties and applications in various areas, including hydrogen storage and utilization in nuclear reactors as neutron moderators [1] and fuel rod cladding materials [2], due to their low thermal neutron absorption cross-section and good mechanical properties. Therefore, numerous experimental and theoretical investigations on their structural and dynamical properties have been conducted and reported in publications since the 1950s.

It is well established that zirconium hydrides have a complex phase diagram as a function of temperature, hydrogen content, and pressure [3,4,5]. Different phases such as ζ-ZrH_0.5 [6], γ-ZrH, δ-ZrH_1.5, and ϵ-ZrH₂ have been identified. A stoichiometric γ-ZrH (face-centered orthorhombic metal sublattice) can be prepared by special temperature treatment [5]. Among those phases, this study focuses on γ-ZrH in the Cccm space group, and ϵ-ZrH₂ in the I4/mmm space group. Though the phase diagram has been determined, some hydrides’ structures, stabilities, or formation mechanisms are still unknown [7].

The optical phonon properties of ZrH_1.5 were first analyzed using inelastic neutron scattering (INS) by Andresen et al. in 1957 [8]. Yamanaka [9] studied the thermal and mechanical properties of zirconium hydrides. Based on first-principles calculations, considerable research efforts have been devoted to the mechanical and thermodynamical properties of zirconium hydrides. Early studies by Ackland found that bistable crystal structures exist with minimal energy difference in ϵ-ZrH₂ [10]. Subsequent studies demonstrated that one of the bistable structures may be unstable or metastable. Elsasser reported that at the center of the Brillouin zone (gamma point), the potential energy profiles for hydrogen atoms in ZrH₂ are parabolic up to about 300 meV [11]. Besides theoretical treatments, INS spectroscopy has also been employed to measure the vibrational spectrum of ϵ-ZrH₂. The peak splitting of multi-phonon events was observed in 1983 by Ikeda [12], yet was not explained. Kolesnikov et al. later attributed those split peaks below multi-phonon events in γ-ZrH and γ-ZrD to bound multi-phonon states [13,14].

Those experiments and calculations revealed that hydrogen vibrations in zirconium hydrides are almost harmonic for the ground and first excited states but turn out to be anharmonic for the higher excitations of the hydrogen atoms in the lattice. Though phonon theory and techniques have been well developed, they are mostly based on the harmonic approximation 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 frequency information can be accessed, which limits further analysis involving the polarization of phonons. More recently, the temperature-dependent effective potential (TDEP) method has been developed to address the anharmonicity problem by taking temperature effects into consideration from MD simulations [15,16]. However, none of the above approaches are capable of obtaining comparable results with our broadband INS data for zirconium hydrides, bringing the necessity of an alternative technique to address the problem.

Based on the above considerations, and to present a combined study regarding the origin of anharmonicity in ZrH_x, we conducted a qualitative and semi-quantitative analysis that accounts for the effects observed in the INS spectra of zirconium hydrides by exploring the potential energy surface (PES) of these materials using density functional theory (DFT) and solving the corresponding eigenfrequencies using quantum theory. The isotope effects were also investigated by checking the γ-ZrD system. In this work, we successfully combined DFT, quantum theory, and INS spectra calculations, to improve our understanding of the mechanical and thermal properties of zirconium hydrides. The methodology can be applied to other metal hydride systems.

2. Methods

Multiple ab initio simulation packages within the framework of DFT have been employed to carry out first-principles calculations in this work, including the Vienna ab initio simulation package (VASP) [17,18], and the Real-space MultiGrid code (RMG) [19,20,21]. 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 [22] and the generalized gradient approximation (GGA) [23] with the parametrization of Perdew–Burke–Ernzerhof (PBE) [24], respectively. A first-order Methfessel–Paxton scheme [25] 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 an ultrasoft pseudopotential and a norm-conserving pseudopotential. The RMG phonon methodology has been recently developed [26], 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 a norm-conserving pseudopotential with exchange-correlation described by the PBE functional was used. Valence electrons for Zr were $4d^{2}5s^2$. The convergence criteria for the root-mean-square (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.

The phonon dispersions of ZrH_x (ZrH_x means one of ϵ-ZrH₂, γ-ZrH, and γ-ZrD) were calculated using the finite displacement method (FDM) as implemented in the Phonopy package [27] employing the forces calculated by both VASP and RMG packages. $2\times 2\times 2$ supercells of ZrH_x (containing 96 atoms for ϵ-ZrH₂ or 64 for γ-ZrH and γ-ZrD) were used with small atomic displacements of 0.01 Å. A gamma-centered $5\times 5\times 5$ k-mesh for ϵ-ZrH₂, and an $11\times 11\times 11$ mesh for γ-ZrH and γ-ZrD were used to integrate total energy. K-mesh and energy cutoff have been tested to converge to 1 meV/atom in RMG and VASP for the total energy. The phonon properties computed by the Phonopy package were used to simulate the INS spectra [28] by the OCLIMAX software [29,30], a program aiming at simulating INS spectra based on vibrational normal modes. A $31\times 31\times 31$ uniformly sampled k-mesh in the BZ was used to represent vibrational modes for our OCLIMAX calculations.

The frozen phonon calculations and potential surface sampling calculations were done with the VASP package.

To validate simulation results, the experimental INS spectrum for ϵ-ZrH₂ was measured at 5 K with the VISION spectrometer [31] at the Spallation Neutron Source (SNS), Oak Ridge National Laboratory (ORNL). The ZrH₂ sample (99% purity) was obtained from Sigma-Aldrich, and for the INS experiment, the sample was placed in a flat aluminum container of $50\times 50\times 0.2$ mm${}^3$ size. The estimated transmission of thermal neutrons through the sample is ∼90%, that means the multiple neutron scattering should be small [32]. The spectrum from an empty container was measured under similar conditions and subtracted from the sample data. INS spectra for γ-ZrH and γ-ZrD were measured previously [13,14] with the TFXA spectrometer [33] at ISIS, Rutherford Appleton Laboratory, UK. INS can measure the vibrational spectrum weighted by the amplitude of motion and the neutron cross-section of the atoms. INS has no selection rules, and more importantly, the comparison of INS simulations with experimental data is rigorous and straightforward. Thus, INS becomes a perfect technique to be combined with theoretical vibrational analysis.

3. Results and Discussion

3.1. Atomic Structures and Phonon Properties of ZrH_x

The crystal unit cells of ZrH_x and their calculated phonon bands are shown in Figure 1. The lattice vectors (a, b, c) used in our simulations for ϵ-ZrH₂ [34] and γ-ZrH(γ-ZrD) [13] were (4.97728, 4.97728, 4.449) and (4.549, 4.618, 4.965) in Å, respectively. Since phonon dispersions (especially optical modes) are sensitive to lattice parameters, no geometry optimization was performed on the lattice vectors of these three materials to better compare simulations with experiments. As is seen in Figure 1, large phonon band gaps exist between high energy modes and low energy modes for all materials.

Figure 2 plots the experimental INS spectra for ZrH_x, and Figure 3 shows the simulated INS spectra for ϵ-ZrH₂ with RMG and VASP along with the VISION data. Within the frequency region below 1000 meV, multiple peaks representing phonon excitations are observed in all materials. Focusing on ϵ-ZrH₂, the peaks with energies lower than 30 meV correspond to lattice modes of heavy Zr atoms in the crystals. The peaks at around 150 meV correspond to single phonon excitations of H atoms from the ground state to their first excited states (fundamental excitations). All higher energy peaks are overtones and combinations, corresponding to multi-phonon events from the ground state to the excited states of higher orders than the first.

The INS spectra simulated from VASP and RMG for ϵ-ZrH₂ show excellent agreement with the experimental data below 200 meV (see Figure 3). However, discrepancies can be observed between the experimental data and the simulated INS spectra in the higher energy range. In the region of overtone peaks, features at around 260, 390, and 510 meV exist in the experimental data, and cannot be observed in the calculation. This behavior suggests that the harmonic approximation, which in this case is good for the fundamental peaks, becomes inaccurate for higher energy excitations; thus, anharmonic effects need to be taken into consideration. It should be noted that for γ-ZrH and γ-ZrD, similar split peaks were also observed in the experimental data (see Figure 2).

3.2. Harmonic Approximation: Frozen Phonon Method

Frozen phonon calculations [38] for phonon modes at the gamma point were performed to explore the hydrogen potential energy profiles in ϵ-ZrH₂ as it usually captures the anharmonicity (if any). The potential energy profiles were sampled along the phonon polarization vectors, and the RMS displacements of phonons ($<u^2{>}^{0.5}$) were chosen to be up to 0.58 Å (corresponding to deformation potential energy at ∼1 eV). Low-frequency lattice modes (smaller than 30 meV) were not included in our frozen phonon calculations.

Our results show that the sampled potentials can be well fitted to parabolas, leading us to conclude that no anharmonic effects can be probed with the frozen phonon method, presumably due to the anharmonicity falls outside the scope of potential energy profiles in polarization directions at the gamma point. To gain more insight into the origin of anharmonicity, we have conducted potential energy surface calculations as presented below.

3.3. A Direct Method: Mapping Eigenfrequencies from Schrödinger Equations with the Simulated INS Spectra

To obtain accurate phonon frequencies, ideally one should treat the collective motion of nuclei in a full quantum approach, i.e., solving the many-body nuclear Schrödinger equation for vibrational excitations. However, this approach is computationally prohibitive: the Schrödinger equation may be easily solved for a single particle, but becomes impossible to solve as the system size significantly increases. In this case, we made three approximations to simplify the complex many-body problem to a single-particle problem.

First, we separate vibrations of H(D) atoms from Zr atoms. The justifications of this separation come from (1) PDOS comparison between γ-ZrH and γ-ZrD with scaled high-energy modes as shown in Figure 4. On the basis that the low-energy modes of γ-ZrD and γ-ZrH are almost identical (the theoretical ratio is 1.005 if Zr and H(D) move in-phase), while the high-energy modes differ by a scaling factor of $\sqrt{2}$ (which is close to the frequency ratio ${\omega }_H/{\omega }_D=\sqrt{m_D/m_H}=1.4136$ if we assume the H(D) and the Zr are almost independent harmonic oscillators), we conclude that Zr and H atoms are almost independent oscillators. This is consistent with the conclusion drawn by Olsson et al. in ϵ-ZrH₂ and γ-ZrD₂ [39]; (2) the charge densities of Zr and H in ϵ-ZrH₂ have near-spherical distribution with slight deformations, which indicates the Zr-H bonds and their couplings are weak [40]; (3) the mean square displacement (MSD) of Zr is much smaller than H(D)’s, as calculated by [28] at 0 K:

$${\mathit{u}}_l^2=\frac{1}{3}·\frac{2.09}{m_l}\sum _q{w}_q\sum _{\nu }\frac{|{\mathbf{ϵ}}_l{(q,\nu )|}^2}{{\omega }_{q,\nu }}$$

(1)

where l is the index of atom, $\nu$ is the index of phonon mode, q is the index of point in BZ, $m_l$ is lth atom’s mass in atomic units, ${\omega }_{q,\nu }$ is the vibrational frequency in the unit of meV, ${ϵ}_l(q,\nu )$ is lth atom’s polarization vector of the $\nu$th mode at point q, and $w_q$ is the normalized weight at point q. The prefactor of 1/3 arises from the powder sample averaging [28]. The MSD was calculated within the low-energy range (<60 meV, excluding the lowest three acoustic modes) and the high-energy range (>60 meV). As can be seen from the MSD calculation results listed in

Table 1. The calculated MSD and $\sqrt{MSD}$ (in brackets) of Zr at (0, 0, 0) and H(D) at $(3/4,1/4,3/4)$ for ZrH${}_x$ in different energy ranges. Units are in Å${}^2$ for MSD and Å for $\sqrt{MSD}$. In $\gamma$-ZrH($\gamma$-ZrD), MSD(Zr) varies slightly (within 3%) between different Zr positions.

Property	<60 meV	>60 meV
$ϵ$-ZrH${}_2$
MSD(H)	0.00063 (0.025)	0.0144 (0.12)
MSD(Zr)	0.0014 (0.037)	0.00000132 (0.00115)
$\frac{MSD\left(H\right)}{MSD\left(Zr\right)}$	0.45 (0.67)	10,909.09 (104.45)
$\gamma$-ZrH
MSD(H)	0.0008 (0.0283)	0.0136 (0.1166)
MSD(Zr)	0.0015 (0.0387)	0.00000079 (0.00089)
$\frac{MSD\left(H\right)}{MSD\left(Zr\right)}$	0.53 (0.73)	17,215.19 (131.21)
$\gamma$-ZrD
MSD(D)	0.0008 (0.0283)	0.0095 (0.0976)
MSD(Zr)	0.0015 (0.0387)	0.000002245 (0.001498)
$\frac{MSD\left(D\right)}{MSD\left(Zr\right)}$	0.53 (0.73)	4231.63 (65.05)

, in the low-energy range, Zr and H(D) atoms move with comparable amplitudes. However, in the high-energy range, the larger ratio of $\sqrt{MSD\left(H\right)/MSD\left(Zr\right)}$ indicates that H(D) atoms are moving with larger amplitudes than Zr atoms at the same frequency. Thus, from the above three results, it can be concluded that the vibrations of H(D) atoms can be viewed as fast degrees of freedom that are decoupled from the dynamics of Zr atoms.

Second, to simplify the problem, instead of looking into all H(D) atoms, we focus on a single H(D) atom, of which the Hamiltonian is

$$\mathbf{H}=-\frac{{\hslash }^2}{2m}{\mathbf{\nabla }}^2+V\left(\mathbf{r}\right)$$

(2)

where ${\mathbf{\nabla }}^2$ is the Laplacian operator on the H(D) atom’s wavefunctions, and $V\left(\mathbf{r}\right)$ is the potential energy at position $\mathbf{r}$ (a set of $V\left(\mathbf{r}\right)$ constitutes the potential energy surface, PES). In our calculations, $V\left(\mathbf{r}\right)$ are the DFT sampled potential energies of the H(D) atom with subtraction of its energy at the equilibrium position ${\mathbf{r}}_0:V\left(\mathbf{r}\right)=E_{DFT}\left(\mathbf{r}\right)-E_{DFT}\left({\mathbf{r}}_0\right)$. By checking the values of the MSD in Table 1, it can be seen that the H(D) atom is mostly vibrating within a sphere of radius 0.12 Å, which is much smaller (∼5%) than the nearest H-H distance (2.24 Å). In fact, interactions of the nearest H atoms are much weaker than the H-Zr interactions—the latter dominates the local potential energy profile and thus the anharmonic behavior of H.

Third, for the sake of simplicity and to present a simple yet robust explanation of the anharmonic effects, we use a 2D representation of $V\left(\mathbf{r}\right)$ and thus $\mathbf{H}$ in Equation (2). It is important to highlight that we do not intend to provide an entirely accurate solution to the anharmonic problem in ZrH_x, but instead, we want to demonstrate the nature of anharmonicity of these systems in a qualitative and semi-quantitative way. The exact 3D many-body problem (for all H(D) atoms) is not only computationally prohibitive but also unnecessarily complicated; however, with the single-particle approximation in 2D, one can easily explore the anharmonicity by visualizing the potential energy surface, and more importantly, solve the single-particle Schrödinger equation and directly compare the result with experiments. Moreover, considering the orthogonality of basis vectors and small displacements of the H(D) atom, it is reasonable to approximate the 3D problem with 2D representations. Even though one may expect a qualitative explanation at best taking the above approximations into account, surprisingly, the results we obtained can semi-quantitatively describe the experiments, and thus validate the approximations made here, as will be further discussed below.

After checking the structure of ϵ-ZrH₂, we decided to sample the (112) plane (shifted so that the equilibrium position of the H atom was included in the plane) as this is one of the planes with the highest symmetry, and involves the most representative shapes of hydrogen PES in ϵ-ZrH₂ in 2D (it contains harmonic potentials in some directions and highly anharmonic potentials in others). The plane was sampled equidistantly with a mesh grid spacing around 0.005 Å (100 points per direction, 10,000 mesh grids totally). Based on the sampled grid, a finer grid with a dimension of $1200\times 1600$ = 1,920,000 was generated using a linear interpolation method to build the Hamiltonian matrix according to Equation (2). A similar scheme has been applied to γ-ZrH and γ-ZrD. Figure 5 presents a view of the plane and the sampled PES.

By numerically solving the corresponding single-particle Schrödinger equation of the H(D) atom, we obtained the eigenfrequencies and wavefunctions, which represent the energy states of the atom, that can be displayed as in Figure 6. The eigenfrequencies of three materials have also been plotted in Figure S2 [41] in the Supplementary Material, and comparisons between eigenfrequencies and experimental spectra can be found in Figure 7. It is particularly important to point out that the eigenfrequencies here are indeed the numerically solved eigenfrequencies subtracted by the ground state frequency, to be comparable with the experimental INS spectra, which are measuring the excitations of phonons from the ground state to high energy states. It can be seen that the anharmonic peak positions in the experimental spectra can be reasonably described by eigenfrequencies for the first three overtones (higher-order overtones in experiments are indistinguishable). For example, in the case of ϵ-ZrH₂ in Figure 7, the eigenfrequencies in red lines at around 264, 384, and 493 meV are below the harmonic multi-phonon bands in magenta lines, and thus correspond to the anharmonic peaks in experimental spectra in black lines. The eigenfrequency at 264 meV agrees well with the experimental value, and the eigenfrequencies at 384 and 493 meV have around 6 and 17 meV differences with the experimental values, respectively. In the cases of γ-ZrH and γ-ZrD, conclusions can be made that the agreements at low energies are pretty good (except for the ∼10 meV difference at around 280 meV in γ-ZrH). It must also be mentioned that in Figure 7c, the ZrD sample was contaminated with H, and thus introduced the extra low-frequency peaks, e.g., at around 140 meV.

The eigenfrequencies from the Schrödinger equation, which represent phonon properties at the gamma point in the reciprocal space, were then convoluted with the scaled energy resolution function of neutron spectrometers in order to better compare the calculations with the experimental spectrum:

$$S_{convoluted}\left(E\right)=\delta \left(E_{eigen}\right)\times \sigma \left(E\right)$$

(3)

where $\sigma \left(E\right)$ is the spectrometer resolution function, E is the energy transfer, $E_{eigen}$ are eigenfrequencies, and $S_{convoluted}\left(E\right)$ is the final convoluted spectrum incorporating anharmonicity from PES. For the indirect geometry instruments (VISION and TFXA) used in our experiments, the energy-dependent resolution function can be approximated by $\sigma \left(E\right)=0.3099+6.1988\times {10}^{-4}E+1.2398\times {10}^{-8}{E}^2$, and both $\sigma \left(E\right)$ and E are in the unit of meV. The convoluted results of ZrH_x are shown in Figure 8.

It can be seen that the convoluted spectra for all materials show reasonable agreement with the experimental INS spectra. Specifically, we can see significant improvements in the prediction of anharmonic peak positions, as both the experimental and the convoluted spectra show similar peaks splitting below the multi-phonon bands. Although there are still some noticeable discrepancies, such as different broadening and intensities of the convoluted spectra compared to the experimental spectra, we think this accuracy is a considerable achievement, given the fact that three prior approximations were made. Moreover, we expect that the broadening and intensity discrepancies can be systematically decreased if we incorporate eigenfrequencies at other k-points than the gamma point into calculations, which, however, is beyond the scope of this work.

To mimic the contamination of the γ-ZrD sample with H in the experiments, H atoms of concentration 3.125 at.% (1 in 32 D atoms in the simulation box was substituted by an H atom) were also added to the pure γ-ZrD in our simulations. The calculation successfully 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). Note that the relative higher intensity of low-energy modes in experimental γ-ZrD rather than simulated low-energy modes of γ-ZrD_0.96875H_0.03125 is mostly due to contributed low-energy modes by α-Zr [13].

Note that the exact knowledge of vibrational dynamics of different Zr-H materials is important not only from the fundamental physics and chemistry points of view, but also for practical applications in the nuclear field such as critical safety studies. Therefore the obtained results on the dynamics of γ-ZrH_x materials in the current study, which describe the anharmonicity of the INS spectra at energies above 200 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{}^{\circ }$C (at which the first overtone level in Zr-H is ∼10% populated).

4. Conclusions

The anharmonicity phenomena in ϵ-ZrH₂, γ-ZrH, and γ-ZrD have been identified and thoroughly studied with DFT and INS techniques. 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, γ-ZrH_x has been considered as a neutron moderator material for nuclear reactors; in this case, the anharmonic effects cannot be ignored since neutrons with high energies will scatter from the hydrogen atoms and will experience deviations from the scattering kernel predictions based on the harmonic approximation. While the harmonic model failed to explain the split peaks shown in the experimental INS spectra, results of eigenfrequencies of the 2D Schrödinger equation convoluted with instrument resolution functions give a good description of the anharmonicity in these materials. One can see obvious anharmonic peaks beyond the harmonic overtones in the spectra, which are in qualitative agreement with the experimental data. The method proposed here, by combining DFT and solutions of the Schrödinger equation compared to INS, is a possible way to explain anharmonicity in materials beyond the harmonic regime. This work could shed light on further vibrational studies in anharmonic systems like TiH_x, PdH_x, and other metal-hydrogen systems.