Mechanical properties of metal dihydrides

The dihydrides considered in this study are all metals, observed in the current calculations, and illustrated in the band structures presented in many previous computational studies of metal dihydrides[13–15, 22]. In metallic systems, the calculations of elastic constants can be profoundly sensitive to the Brillouin Zone sampling and the effective electronic temperature [23, 24]. In bcc-Ta, the equilibrium lattice parameter and total energy converged relatively quickly with k-point density, but the elastic constants, and particularly the shear moduli, C₄₄ and C_s, defined as

$$\begin{eqnarray}&&{{C}_{s}}=\frac{\left({{C}_{11}}-{{C}_{12}}\right)}{2}\end{eqnarray} \tag{ 1 }$$

showed large dependence with k-sampling (5 GPa variability) out to k-points samples as large as 40³ [24]. Hence, while converged calculations of the ground state lattice constant, total energy, and bulk modulus are relatively straightforward to obtain, computation of elastic constants might have much greater uncertainties and need to be approached more skeptically.

Figure 2 illustrates the k-sampling dependence of the bulk modulus and the elastic constants computed for cubic ErH₂, from a 8³ through a 24³ k-sampling of the full Brillouin Zone for the 3-atom primitive cell. In all cases we consider, the lattice constant and total energy exhibit almost no k-point dependence, i.e. they have converged to within 0.001 Å and 2 meV/formula-unit (f.u.) over the range of k-point samples investigated. The bulk modulus also converges very well (to within 0.1–0.2 GPa) and relatively quickly. The other elastic constants exhibit much greater variability, as seen in earlier analyses of elastic constant computations in metals [23, 24]. Once the sampling gets denser than 12³, the computed constants appear to converge around an average value, and there is little apparent diminution in the scatter of the computed values about this average out to the largest k-point sampling we consider. This variability with k-sampling dominates the uncertainties associated with computation of elastic constants and is the principal factor in the numerical precision with which the elastic constants can be predicted. The typical numerical uncertainties due to k-sampling for computation of shear moduli in the cubic dihydrides are 2–3 GPa, but are as large as 5–6 GPa in some cases.

Zoom In Zoom Out Reset image size

It is pointless to attempt to numerically refine these elastic constants any further. First, the apparent rate of convergence, narrowing of the width of the scatter, is very slow with number of k-points but, more importantly, this k-dependent scatter has become smaller than the expected accuracy of the DFT. The B computed with semi-local functionals such as PBE can typically only predict experimental results for a bulk modulus to within 10–20%. The rapidly increasing computational cost of using denser k-sampling does not improve the total quality of the numerical predictions, i.e. the increased computational cost is not recouped in increased fidelity when these physical uncertainties are also considered.

This variability does not indict the quality of the calculation of C₁₁ and C₁₂ at each k-sampling density (from which the elastic modulus C_s is constructed). The bulk modulus in a cubic system, is also constructed from these same elastic constants:

$$\begin{eqnarray}&&B\left[\text{cubic}\right]=\frac{{{C}_{11}}+2{{C}_{12}}}{3}.\end{eqnarray} \tag{ 2 }$$

The stress-derived B({C_ij}) show very small variability (typically  ∼0.1 GPa) with k-point sampling. As seen in figure 2, C₁₁ and C₁₂ suffer significant variations with k-point sample out to the largest sampling we consider, but these values are very closely correlated, leading to an almost constant bulk modulus. An independent calculation of the bulk modulus, B(E) from a fit to a third-order Birch–Murnaghan equation of state of the total energies at different lattice parameters, yields the same result for the bulk modulus as the stress-derived calculation, to within  ∼1 GPa, verifying this evaluation of B({C_ij}) from the elastic constants. Moreover, the lattice constant and B from the stress-derived calculation and from the total energy Birch–Murnaghan equation of state closely match, indicating the (artificial) electronic temperature used to facilitate self-consistency in these metallic calculations gives a close approximation to a 0 K limit. This isolates the k-dependence as responsible for the variation in evaluated C_ij, and not numerical artifacts from other aspects of the computational model.

Figure 3 illustrates the convergence of the elastic constants with k-sampling for the tetragonal structure HfH₂. In tetragonal dihydrides, the k-dependence of bulk properties becomes more severe, and spreads to the assessment of the lattice parameters. This can be attributed to a sensitivity of the structure and energy to occupation of states at the Fermi level. The associated computational challenge is to resolve the larger DOS at the Fermi level in the tetragonal dihydrides [15]. As in the cubic structure, the shear elastic constants are particularly sensitive. While the lattice constants a₀ and c₀ (not plotted) show greater variability, the equilibrium volume converges better with k-point sampling, converging as well as the volume (and lattice constant) had for the cubic crystal.

Zoom In Zoom Out Reset image size

The bulk modulus can be evaluated in terms of the elastic constants as

$$\begin{eqnarray}&&B\left[\text{tetragonal}\right]=\frac{\left({{C}_{11}}+{{C}_{12}}\right){{C}_{33}}-2C_{13}^{2}}{{{C}_{11}}+{{C}_{12}}+2{{C}_{33}}-4{{C}_{13}}}\end{eqnarray} \tag{ 3 }$$

in a tetragonal crystal [39]. This stress-derived B({C_ij}) also shows more rapid convergence with k-sampling, although not quite as fast nor quite as tight as for the cubic structure.

The bulk modulus can alternatively be computed from the equation of state. The B(E) can be obtained from analytic derivatives of a Birch–Murnaghan fit to the total energies of a sequence of pressure-optimized tetragonal cells. The energy-fit B(E) match very well the stress-derived B({C_ij}), to within the same  ∼1 GPa that was observed for the cubic crystal, verifying that the calculation of elastic constants using strain-derivatives of the stress is consistent with the total energy calculations for the tetragonal structure. This verifies the form and evaluation of equation (3), and verifies the computation of the elastic constants within the formula given in equation (3).

The results for the structure and elastic constants for the cubic structure of ErH₂ and the Group-IIIB and Group-IVB dihydrides are presented in table 1. The k-limit values are the full converged prediction of the DFT simulations. The quoted uncertainties in the k-limit values reflect numerical uncertainty in the converged average value, taken as the magnitude of the scatter from the average as the k-point sampling is increased. The value using a 12³ k-sampling is also presented, a representative of the value one obtains with just a single k-point sample, at the threshold at which the simulations appear to converge to a k-limit asymptotic average. The 12³ k-sample results in table 1 are at the boundary of the asymptotic range, i.e. converged within the uncertainties of the k-limit values. The formation energy and lattice constants are fully k-sampling converged, to the precision quoted in table 1.

Table 1. Structural properties and elastic constants computed for cubic metal dihydrides.

	$-\Delta {{E}_{f}}$ (eV)	a0 (Å)	C11 (GPa)	C12 (GPa)	C44 (GPa)	Cs (GPa)	B({Cij}) (GPa)	B(E) (GPa)
ErH2	2.281a	5.126b
k-limit	—	a0(expt.)	137(5)	61(3)	76.8(0.5)	38(5)	86.5(.1)	—
k-limit	2.071	5.113	140(4)	63(2)	77.8(0.5)	39(3)	89.1(.1)	88.3(.1)
k  =  123	2.075	5.113	145.4	61.1	78.5	42.2	89.2	88.4
ScH2	2.073a	4.783c
k-limit	1.990	4.766	166(5)	63(2)	80(1)	51(4)	97.4(.2)	97.3(.1)
k  =  123	1.992	4.766	168.8	61.7	76.9	53.6	97.5	97.4
YH2	2.346a	5.195c
k-limit	2.091	5.206	132(3)	58(2)	73(1)	37(3)	82.8(.1)	82.9(.1)
k  =  123	2.095	5.205	132.6	56.2	73.6	38.2	82.8	82.9
LaH2	2.151a	5.663d
k-limit	1.697	5.710	95(2)	48(1)	53(1)	23(1)	63.5(.1)	63.8(.1)
k  =  123	1.698	5.710	95.7	47.5	54.1	24.1	63.6	63.9
TiH2	1.283a	4.475e,4.463f
k-limit	1.434	4.417	119(1)	155(1)	−21(4)	−18(1)	143.1(.1)	140.2(.1)
k  =  123	1.435	4.417	120.5	154.3	−21	−17	143.0	140.3
ZrH2	1.687a	4.794g
k-limit	1.614	4.811	88(2)	157(1)	−30(4)	−34(2)	133.9(.1)	134.4(.1)
k  =  123	1.615	4.810	89.0	156.3	−30	−33.7	133.9	134.5
HfH2	1.375a	4.713h
k-limit	1.376	4.713	88(3)	179(2)	−42(3)	−45(3)	148.5(.1)	149.7(.1)
k  = 123	1.376	4.713	81.5	182.1	−49	−50.3	148.6	149.6

^aExperimental result from Beavis [40]. ^b[41]. ^c[42]. ^dAt 0 K, [43]. ^eAt 300 K [44] (N.B., ground state is tetragonal below 300 K [45]). ^fExtrapolated to TiT₂, from sub-stoichiometric data using linear relationship in [46]. ^gExtrapolated to ZrH₂ from sub-stoichiometric data using linear relationship in [47]. ^hFor HfH_1.7, from [48]. Note: The value in the parentheses for the k-limit values reflects its k-sampling uncertainty, variation with respect to k-sampling in the large k-sampling limit.

The elastic properties are not much changed for ErH₂ if computed at the experimental lattice parameter rather than the computed a₀. Usually the PBE functional significantly overestimates the lattice constant and therefore strongly underestimates the resulting computed bulk modulus, and elastic properties computed at the experimental lattice parameter improve the accuracy of the calculations. For the rare earth dihydrides, the bare PBE lattice constant, absent zero-point effects and magnetism, yields exceptionally good agreement with experiment, PBE only  <0.3% different and smaller. This fortuitously good agreement in a₀ means the predictions lack a bias from a mistaken volume. Nonetheless, that the small 0.3% change in a₀ leads to almost a 3 GPa change in the computed bulk modulus serves as a cautionary result, slight changes in the structure could result in significant changes in dependent properties such as the bulk modulus, injecting another source of error into the analysis.

The formation energy $\Delta {{E}_{f}}$ for a dihydride of metal $\text{M}$ is obtained from the computed total energy of $\text{M}{{\text{H}}_{2}}$, and of the ground state of the bulk metal and hydrogen molecule, ${{E}_{\text{M}}}$ and ${{E}_{{{\text{H}}_{2}}}}$:

$$\begin{eqnarray}&&\Delta {{E}_{f}}={{E}_{\text{M}{{\text{H}}_{2}}}}-{{E}_{\text{M}}}-{{E}_{{{\text{H}}_{2}}}}.\end{eqnarray} \tag{ 4 }$$

Thermal effects are neglected in equation (4), they are minor on the scale of errors in the physical approximation of the DFT. Zero-point effects are also ignored. The largest zero-point correction to the formation energy reduces the formation energy by the zero-point vibration of the hydrogen molecule [49]. This molecular constant is partially counterbalanced by the zero-point vibrational energy of the hydrogen atom in each metal lattice. The H₂ ZPE gives a sense of the magnitude of the remaining errors in E_f and could be trivially added to the computed values for improved comparison to experiment. The ZPE had been computed in the PBE simulation context to be 0.224 eV [49] for this ZPE, but we note that it would be a more accurate and appropriate physical comparison to use the 'exact' ZPE for H₂, 0.270 eV/H₂ [50] in any such adjustment. The agreement of our computed formation energy, ΔE_f in table 1, with experimental heats of formation, as given by Beavis [40], is surprisingly good for ErH₂ and ScH₂ and still reasonable for the other hydrides, despite neglecting thermal and quantum nuclear effects. The ΔE_f calculated here for ErH₂ matches a VASP-calculated value to within 3 meV [49].

The computed lattice constant for β-ErH₂, 5.113 Å, agrees well with experiment [41], 5.126 Å. Our VASP calculation gives a₀  =  5.129 Å, consistent with this result. This agreement between the two codes, with very different pseudopotentials and basis sets, provides additional verification of the computational methods. Our results also agree with previous VASP calculations [49], who quote a₀  =  5.128 Å. The computed lattice parameters for the other dihydrides presented table 1 agree well with experiment [42–48] and with previous calculations [18, 21, 22, 51–55]. The stoichiometric group-IVB dihydrides are known to be tetragonal, and adopt the cubic phase for a small range of substoichiometric H. The TiH₂ transforms into a fcc-CaF₂ phase above 310 K [45], but adopts the tetragonal structure at lower temperatures. In table 1, the first quoted experimental lattice parameter for TiH₂ is for this observed high-temperature cubic phase [44]. The substoichiometric hydrides show only a small dependence of lattice parameter on the H content over the range of H content where the cubic phase is stable, and the relationship between H content and lattice constant is roughly linear. The second quoted 'experimental' lattice parameter for TiT₂ is obtained by extrapolating to the stoichiometric ditritide using the relationship developed [46] over the range of substoichiometric tritide where the fcc phase is stable (1.5  <  T/Ti  <  1.8). This extrapolated a₀ in TiT₂ is slightly smaller, partially attributable to an isotope effect, the hydrogen vibrations causing a larger expansion in the lattice than the tritium.

Unlike for TiH₂, no cubic form of the stoichiometric ZrH₂ has been observed, but, similarly, a linear relationship was developed for the substoichiometric ZrH_2−x system [47], which we use to extrapolate to obtain the 'experimental' value of a₀  =  4.794 Å quoted in table 1 for the dihydride. The lattice constant measured for HfH_1.7 is used as the experimental value here for the dihydride; we were unable to find the experimental data to perform a comparable extrapolation from substoichiometric data for HfH₂. The agreement of the calculated a₀ with these inferred 'experimental' lattice constants for the cubic structure of group-IVB dihydrides is very good.

The C₁₁, C₁₂, and C₄₄ elastic constants are computed from the strain-derivatives of the stress, the shear modulus C_s and B({C_ij}) are evaluated from these using the analytical relationships in equations (1) and (2). The alternate energy-fit calculation of B(E) agrees with B({C_ij}) across this set of dihydrides, verifying the computation of B. The computed bulk modulus for TiH₂ is less good, the difference is slightly larger, 3 GPa. This signals a potential numerical issue in the calculation, and upon investigation we discovered that this larger discrepancy was due to inaccuracy in the real space integration in the total energy evaluation. A finer real space grid, increased from 24³ to 32³ (or more) brought the B(E) verification into the same agreement as for the other dihydrides, and had minimal effect on the elastic constants computed from the stress-strain relations. This exercise illustrates how seemingly redundant verification checks, with meaningful estimates of expected and actual uncertainties, can identify potential inaccuracies in the calculation, and thereby add greater confidence to the computational results.

Quijano, et al [22] had done a notably careful study of the structural properties of the group-IVB dihydrides, examining stability of the fcc versus the fct phase, additionally examining the importance of relativistic effects. Their results for fcc phase TiH₂, ZrH₂, and HfH₂ lattice parameters of cubic (4.428, 4.817, and 4.727 Å, respectively) and the bulk modulus (144, 136, 148 GPa) are closely mimicked in the current results. We find that these group-IVB cubic dihydrides all have negative shear moduli, both the C₄₄ and the tetragonal shear C_s, indicating that the cubic structure is unstable to a distortion. Given that these are all known to take tetragonal ground states, this is not surprising, but the result indicates that the cubic phase does not even have a locally stable region, and, further, suggests that these shear elastic constants might be used computationally to detect the instability of the cubic phase, a notion we examine later. This structural instability leads to numerical instabilities in the assessment of the stresses in the calculation of the C₄₄. The computed stresses are negative, and not cleanly linear, so the numerical evaluations of C₄₄ from the stress-strain results are much more uncertain at each k-point.

We now shift to erbium's immediate neighbors, to examine the mechanical properties of cubic dihydrides in the rare earths (RE). A search was performed on the Pearson database [56] of crystals to identify all the RE dihydrides that formed a cubic (Fm$\bar{3}$m) dihydride ground state. The cubic structure is the stable ground state across almost the entire lanthanide series. The ground state crystal structure of PmH₂ is not known. Promethium has no stable isotope and only trace amounts are found in naturally occurring ores (it has never been made). The dihydrides of Eu and Yb instead crystallize into an orthorhombic structure (Pnma, #62), rather than the cubic structure of ErH₂ and the other RE dihydrides. This could be anticipated, as Eu and Yb do not always exhibit the same [+3] oxidation state that characterizes the other rare earths, but instead can also adopt a [+2] oxidation state, i.e. with respect to the other RE atoms, they will have one fewer valence electron.

Table 2 presents the DFT results for the lattice parameter and elastic constants computed for the RE dihydrides in the fcc CaF₂ structure. Results are presented both from SeqQuest calculations using a linear combination of atomic orbitals (lcao) basis set (using the 3-atom primitive cell and a 20³ k-point sampling), and also from calculations using VASP with a plane wave (pw) basis. In these results, to examine the trends across the RE series, the lcao calculations for EuH₂ and YbH₂ use pseudpotentials with a [+3] oxidation state (rather than the [+2] oxidation state pseudpotential), and in the cubic (fcc) structure rather than the Pnma orthorhombic structure observed in nature.

Table 2. Lattice parameters a₀ (Å), elastic constants (GPa), and computed Debye temperature (K) for cubic RE dihydrides.

Hydride		a0(Exp.)a	a0(DFT)	C11	C12	C44	B({Cij})	Cs	${{\Theta}_{\text{D}}}$ (DFT)d
LaH2	lcao	5.667	5.710	93.5	48.5	53.7	63.5	22.5	364
	pw		5.658	89.3	48.5	47.3	62.1	20.4	343
CeH2	lcao	5.581	5.622	98.6	50.6	56.2	66.6	24.0	369
	pw		5.437	90.7	55.7	46.7	67.4	17.5	324
PrH2	lcao	5.518	5.552	102.2	51.6	58.4	68.5	25.3	374
	pw		5.555	102.0	52.1	53.1	68.7	25.0	363
NdH2	lcao	5.464	5.487	106.3	52.8	60.6	70.7	26.7	376
	pw		5.490	107.4	54.0	55.5	71.8	26.7	366
PmH2	lcao	b	5.428	111.8	54.5	63.1	73.6	28.7	382
	pw		5.431	111.6	55.6	58.3	74.3	28.0	372
SmH2	lcao	5.374	5.375	116.8	55.5	65.5	75.9	30.6	383
	pw		5.385	115.8	56.3	60.2	76.1	29.8	372
EuH2	lcao	c	5.325	121.8	56.4	67.8	78.2	32.7	388
GdH2	lcao	5.303	5.279	126.8	57.5	70.4	80.6	34.7	388
	pw		5.285	122.3	58.7	65.4	79.9	31.8	374
TbH2	lcao	5.246	5.234	130.9	58.2	72.4	82.4	36.3	391
	pw		5.237	127.0	59.5	68.1	82.0	33.8	379
DyH2	lcao	5.201	5.191	135.4	59.5	74.3	84.8	37.9	392
	pw		5.200	131.7	62.0	68.7	85.2	34.8	377
HoH2	lcao	5.165	5.151	139.5	60.2	76.8	86.6	39.7	395
	pw		5.160	133.7	62.6	73.5	86.3	35.6	382
ErH2	lcao	5.123	5.113	144.4	61.3	78.2	89.0	41.6	396
	pw		5.122	137.5	63.7	74.6	88.3	36.9	382
TmH2	lcao	5.090	5.074	147.5	62.4	81.2	90.8	42.5	399
	pw		5.090	140.1	66.0	75.4	90.7	37.1	381
YbH2	lcao	c	5.039	150.3	63.3	82.7	92.3	43.5	397
LuH2	lcao	5.033	5.004	153.7	64.5	84.4	94.2	44.6	398
	pw		5.020	145.7	66.6	78.5	93.0	40.0	381

^aExperimental lattice parameter for LaH₂ and CeH₂ from [57]; PrH₂ through LuH₂ from [58], except GdH₂ from [59]. ^bPmH₂ has never been made. ^cComputed result shown is for a pseudopotential with the [3+] oxidation state. ^dComputed from elastic constants using relation from [60].

The agreement between the lcao and pw calculations for each RE dihydride is very good across the series. The differences in elastic constants between the pw and lcao results are consistent with the numerical uncertainties with respect to k-sampling reported in table 1. The bulk modulus, as illustrated in figure 2 for ErH₂ and analyzed in table 1 for several dihydrides, exhibits much finer convergence with k-point sampling, so that the total uncertainty in this quantity is dominated by the form of the equation of state, estimated above to be  ∼1 GPa, rather than k-sampling dependence. The differences between the lcao and pw results for B in the RE series are typically  ∼1 GPa, consistent with the estimated uncertainty in the model form used to extract the elastic properties. Because the k-point sampling used in computing the elastic constants for the RE series are well within the asymptotic range used for computing the k-limit result, the uncertainties (with respect to the k-limit value) due to k-sampling in the values quoted in table 2 will be the same as those presented for the stable cubic dibidrides in table 1, i.e.  ∼4 GPa for C₁₁ and C_s, ∼2 GPa for C₁₂, ∼3 GPa for C₄₄. The values in table 2 are quoted to tenths of GPa only to facilitate finer comparisons between the lcao and pw calculations, and also to analyze trends across the series, this finer precision is otherwise physically meaningless. That the results using very different pseudopotentials and basis sets, using different methods for extracting elastic constants for the results, and using different cells and k-point sampling, are in good agreement lends greater confidence to the numerical predictions. The different approaches give chemically and numerically consistent results.

The computed properties of the RE dihydrides are similar to one another and follow a monotonic trend from one end of the lanthanide series to the other. The small exception is the plane wave result for C₁₂ (and the C_s derived from it) for CeH₂. Ce is another rare earth that often adopts an oxidation state different from [+3]: it can promote its f-electron into the valence electrons to access a [+4] oxidation state. The lcao result uses a pseudopotential leaving exactly one f-electron in the core, imposing the [+3] oxidation state that is consistent with its neighbors, while the PAW pseudopotential of the PW calculation does not make this constraint on the oxidation state.

The lattice parameter steadily shrinks across the series, the effective size of the RE atom getting smaller as the nuclear charge increases. The computed a₀ are in very good agreement with experiment [43, 58, 59]. The PBE calculation gives a₀ slightly larger than experiment early in the series and then trends to become slightly smaller than experiment at the end of the series.

All the computed shear moduli are positive, indicating that the RE dihydrides are (at least locally) stable in the cubic phase, consistent with the experimental observations [56]. The elastic constants all steadily become stiffer, as would be expected as the lattice parameter decreases for a given structure. This supports the notion that the bonding and electronic structure are very similar across the RE dihydrides, the dominant difference being the shrinking in size in the RE atom across the series. Overall, DFT with the PBE functional appears to give a very good description of the mechanical properties of cubic dihydrides. Using the relations developed by Deligoz et al [60], an estimate for the Debye temperature can be obtained directly from the elastic constants. The computed results are presented in the last entry in table 2. Unremarkably, the variation across the lanthanides is small.

The negative shear elastic constants for the group-IVB dihydrides presented in table 1 signal that the cubic phase is unstable to a distortion for these materials. Table 3 summarizes the structural properties of the tetragonal (fct) ground state group-IVB dihydrides. In agreement with experiment[61–63], and recent theory [21, 22], the tetragonal form with c/a  <1 is the ground state for all of these, favored over the c/a  >  1 form. The preference for the c/a  <  1 is weak, only 4 meV f.u.⁻¹ in TiH₂ and 12 meV f.u.⁻¹ in both ZrH₂ and HfH₂, but meaningful within the numerical uncertainties present in these calculations. These lattice parameters and energy preferences match well with the relative structural distortion energies computed by Quijano, et al [22].

Table 3. Structural properties for tetragonal metal dihydrides.

		$-\Delta {{E}_{f}}$ (eV)	a0 (Å)	c0 (Å)	c/a	V(fct) (${\mathring{\text{A}}^{3}}$)	V(fcc) (${\mathring{\text{A}}^{3}}$)	$\Delta$E(fcc) (meV f.u.−1)
TiH2	Exp.a		4.528	4.279	.945
c/a  <  1	k-limit	1.445	4.531(2)	4.184(4)	0.924(1)	21.47	21.54	−11
	k  =  123	1.445	4.524	4.198	0.928	21.48	21.54	−10
c/a  >  1	k-limit	1.441	4.329(1)	4.591(2)	1.061(1)	21.51	21.54	−7
	k  =  123	1.440	4.331	4.587	1.059	21.51	21.54	−5
ZrH2	Exp.b		4.9825(9)	4.4488(9)	0.8929(8)
	Exp.c		4.975(3)	4.447(3)	0.894(1)
c/a  <  1	k-limit	1.646	5.011(3)	4.400(6)	0.878(1)	27.62	27.82	−32
	k  =  123	1.646	5.018	4.385	0.874	27.61	27.83	−31
c/a  >  1	k-limit	1.634	4.647(1)	5.135(3)	1.105(1)	27.73	27.82	−20
	k  =  123	1.634	4.644	5.142	1.107	27.73	27.83	−19
HfH2	Exp.d		4.919	4.363	0.886
c/a  <  1	k-limit	1.418	4.935(3)	4.267(5)	0.865(2)	25.98	26.17	−42
	k  =  123	1.418	4.942	4.254	0.861	25.98	26.17	−42
c/a  >  1	k-limit	1.406	4.525(2)	5.094(3)	1.126(2)	26.08	26.17	−30
	k  =  123	1.405	4.520	5.107	1.130	26.08	26.17	−32

^a[45]. ^bAt 293 K, [61]. ^cAt 294 K, [62]. ^dAt 300 k, [63].

We performed additional calculations in the LDA, and obtained the same structural preference for TiH₂ and ZrH₂, of roughly the same magnitudes as with the PBE functional. The result by Ackland favoring c/a  >  1 for ZrH₂ [19] is probably attributable to an inadequate k-point sampling (the quoted k-sampling [19], a sampling less than the 12³ asymptotic threshold for reliable results determined in this study). The result favoring c/a  >  1 for TiH₂ by Wolf and Herzig [15], using an all-electron method, is more mysterious. But the energy preference in TiH₂ involved is so small, 3–4 meV f.u.⁻¹, and may not be physically significant; this is perhaps within the numerical accuracy that can reasonably be expected of older computational DFT methods.

Negative shear elastic constants signal the instability of the cubic structure to a distortion, but do not predict the form that this distortion will take. Given that both the planar shear and tetragonal shear constant presented in table 1 for the group IVB fcc dihydrides are negative, any candidate distortion is plausible. A negative C_s alone would be suggestive of a tetragonal distortion, but that a negative C₄₄ accompanies it indicates other distortions would also be energy-lowering. We investigated rhombohedral crystal distortions, compressing and expanding the cubic structure along the (1 1 1)-axis, for both TiH₂ and ZrH₂. Both compression and expansion lead to energy-lowering distortions, and these hexagonally strained structures are favored over the fcc structure just as both tetragonally-strained fct structures are favored over the fcc structure.

The ZrH₂ favors the (1 1 1)-compressed rhombohedral structure over the fcc structure by 7 meV, and the (1 1 1)-expanded structure is favored by 5 meV, somewhat smaller than energy-lowering obtained with the tetragonally distorted structures. The TiH₂ also favors the (1 1 1)-compressed hexagonal structure over (1 1 1)-expanded structure, at 6 and 1 meV f.u.⁻¹ below the cubic structure. We note that the compressed hexagonal structure is actually more stable than the c/a  >  1 tetragonal distortion. These results suggest that the cubic fcc structure is a local maximum in the structure, surrounded by a valley of lower energy distortions in all directions. This would more readily facilitate the fct  →  fcc transition observed at 310 K [45] for TiH₂, the barrier to a homogeneous distortion would be lower than a barrier represented by a transition through the cubic structure.

The lattice parameters for the fct structures show slightly more variability than for the cubic structure with k-point sampling. Whereas the cubic structures converged to less than 0.001 Å at a 12³ k-sample, the k-sampling variability in the optimized a₀ and c₀ for the fct structure can be as much as 0.2%, a numerical uncertainty that is as large or larger than the uncertainties in the experimental measurements. However, this variation in the computed lattice parameters a₀ and c₀ is tightly correlated. The total cell volume is converged to better than 0.01 ${\mathring{\text{A}}^{3}}$ already at 12³, indicating that the details of the sampling of the electronic states at the Fermi level have a more significant effect on the shape of the structure, the c/a ratio, than its volume.

Table 4 presents the computed elastic constants of the ground state (c/a  <  1) tetragonal metal dihydrides. The uncertainties in the evaluated elastic constants due to k-sampling are comparable, if slightly larger than for the cubic structure presented in table 1 (leaving aside the pathologies that manifest in the calculation of C₄₄ due to the instability of the cubic structure). The tetragonal crystal has a lowered symmetry from the cubic crystal, there are more independent elastic constants to compute than for the fcc structure, and the calculations are more computationally demanding, making it more challenging to verify the results. A couple of the consistency checks we performed are presented in table 4.

First, we computed C₁₃ and C₃₁ independently. To obtain C₁₁, an axial strain along the c-axis is used, and the C₁₃ can be obtained from the derivative of the off-axis stress within the same strain. To obtain C₃₃, it is necessary to perform a strain along the a-axis. The C₃₁ can be obtained from a stress-derivative along a different axis than this strain. By symmetry, C₁₃ and C₃₁ should be identical, but numerically will differ due to different real-space grid spacing and reciprocal space grid spacing along the direction of the strain. This difference can be used as a consistency check, to assess how well the elastic constant calculation is numerically converged with respect to the integration grids. The difference between C₁₃ and C₃₁ is  ∼1 GPa or less for all the tetragonal dihydride calculations, less than the expected k-sampling uncertainty for any individual elastic constant.

A second consistency check involved independent calculations of the bulk modulus. The bulk modulus for these tetragonal crystals was first obtained from the computed elastic constants using the analytic relationship of equation (3). For the 12³ k-sample, we also obtained the energy-fit bulk modulus B(E) from an analytic derivative of a Birch–Murnaghan equation of state, fit to a series of optimized cells under a series of compressive and tensile pressures. This more holistic test assesses both the consistency of the energy and stress calculations and also the efficacy and proper application of the relationship between the elastic constants and the bulk modulus for a tetragonal crystal given in equation (3). The agreement between the stress-derived and the energy-fit bulk modulus is better than 1 GPa, once more within the expected accuracy of the elastic constant calculations, and also matches the known uncertainty in the numerical fit to the equation of state.

For ZrH₂, there are two earlier sets of results for the elastic constants using PBE [52, 54], that are reasonably well-converged and compatible with the current study. Our computed elastic constants match the previous results of Zhang et al [54], and Olsson, et al [52] reasonably well. Our calculations were derived from stress-strain relationships, and these authors instead used a series of energy-fits to obtain their elastic constants, which might explain the modest differences between the different sets of results. All the IVB dihydrides pass the stability tests for the tetragonal structure, and, despite differences in an individual elastic constant as large as 10 GPa (in C₁₁), the resultant bulk modulus still agrees well between these different DFT calculations.

These metals tend to form hydrides over a wide range of H stoichiometry; for instance, erbium and the other RE elements can be readily hydrided to form up to a trihydride [40]. The stability of the different erbium hydrides, ErH, ErH₂, and ErH₃, was investigated with first-principles calculations [49]. The interest in erbium for tritium storage is in the cubic structure. The single phase cubic form is stable for Er:H ratios ranging from 1.9 to 2.1 [40]. In the tritide, the stoichiometry of the initial manufactured sample will change with aging, as tritium decays into helium and modify the mechanical properties of the material. Helium replaces hydrogen in the lattice, and then the helium diffuses away and coalesces into bubbles, leaving behind a substoichiometric tritide. These compositional changes can potentially affect the elastic properties (and the stability of the cubic phase) and therefore are important for modeling the mechanical behavior of aging tritides.

To examine the trends of the elastic properties with a change in stoichiometry, we adopt an averaged-atom approach, where every hydrogen H in the lattice is replaced with a hydrogen atom H^z with a modified nuclear charge z. The advantage of this averaged-atom approach is that it preserves the full symmetry of the crystal in the small primitive cell while exploring off-stoichiometry effects. For $z~<~1$, this approach mimics the average population of x  =  2(1  −  z) vacancies on the H site in the structure, i.e. a random ErH_2−x structure, while for $z~>~1$, this approach mimics a random ErH_2−xHe_x structure where x  =  2(z  −  1). The results for the elastic properties of ErH^z are presented in table 5 for z  =  0.75 to z  =  1.50, using the theoretical lattice constant appropriate to each H^z.

Table 5. Computed structural and elastic properties (in GPa) for cubic ErH$_{2}^{z}$ (k  =  20³).

	a0(Å)	C11	C12	C44	Cs	B({Cij})
ErH$_{2}^{0.75}$	5.226	124.1	49.0	58.8	37.5	74.0
ErH$_{2}^{1.00}$	5.113	144.4	61.3	78.2	41.6	89.0
ErH$_{2}^{1.25}$	5.060	117.9	80.3	73.1	18.8	92.9
ErH$_{2}^{1.50}$	5.088	82.6	105.5	−39.1	−34.3	82.6
ErHHe	5.238

The averaged-atom results indicate that C₁₁ decreases and C₁₂ increases with replacement of H with He in the lattice, culminating with ${{C}_{11}}<{{C}_{12}}$ for ErH$_{2}^{1.5}$, indicating a structural instability. The shear moduli C₄₄ and C_s decrease and then sharply drop at z  =  1.5. The optimized tetragonal structure for ErH$_{2}^{1.5}$ is more stable by 57 meV f.u.⁻¹ than the cubic crystal. Given that ErH^1.5 is isovalent—the same number of valence electrons—with HfH₂, this is not surprising. The lattice parameter first decreases and then turns slightly larger with increasing He concentration. The bulk modulus does not change much across this composition change; only the shear moduli exhibit strong changes.

The numerical results for the averaged-atom structures should not be taken too literally. An ErHHe fcc structure (where every second H atom in ErH₂ is replaced by a He atom), isovalent with the ErH$_{2}^{1.5}$, has a significantly larger lattice constant than ErH$_{2}^{1.5}$. This cubic ErHHe model, nonetheless, exhibits the same structural instability as ErH$_{2}^{1.5}$, although with a larger margin, 157 meV, favoring the tetragonal distortion.

The effects of transmutation of H into He in the lattice have been studied previously [64] using an explicit-atom approach. Using a larger supercell (96 atoms), between one and six H atoms in a larger 96 atom ErH₂ supercell were replaced with He atoms, probing a gradation in He composition much finer than considered here. For small x, those results show C₁₂ decreases rather than increases, and the bulk modulus decreases significantly, and over a much smaller gradation of He composition. This indicates that effects for small concentrations of clustered He on H lattice sites the mechanical properties will be different than for dispersed concentrations of He.

Results for dispersed He should be taken as illustrative of trends rather than as physically predictive. The tritides are not observed to sustain even small concentrations of He located on H sites, but rather the He generated from decayed tritium diffuses away and segregates into bubbles almost immediately [9]. The macroscopically observed changes in mechanical properties of the aging material are dominated by microstructural effects: the bubble density, shape, and size, the pressures within the He bubble interacting with the elastic properties of the remaining tritide matrix, and the dislocations and mesoscale features this interaction generates [3]. The components that go into such an analysis include the elastic properties of the matrix, but that matrix will retain minimal residual He concentrations. The computational results are more useful, indicative of incipient local structural instabilities, if a few He atoms, before they diffuse into bubbles, are found within close proximity of each other. This instability to local distortions due to heterogeneity of local He concentrations could perhaps facilitate increased diffusion of He, and, more speculatively, the statistical likelihood of a near approach of several He might be a contributor to the initial nucleation of bubbles. Tuning ErH$_{2}^{z}$ to z  =  1.5 creates a system isovalent with HfH₂. If several H are replaced by He in a local region, this might favor local distortions and strains.

Simulations of the substoichiometric tritide, on the other hand, should provide more realistic models of reality, either for a tritide that was less than fully loaded, or because He segregation in the aging tritide leaves behind a substoichiometric matrix. The results in table 5 indicate that with an average of 1/4 H vacancies, the lattice becomes slightly (2%) larger and all the elastic constants get slightly softer (∼20%)—relatively small changes considering the large composition change. XRD studies of aging ErT₂ show a linear anisotropic increase in lattice constant of  ∼1% when 1/8³H have converted to ³He [65], consistent with this result. For use in a macroscopic model of the mechanical response of the aging material, this kind of variability in the elastic properties in the matrix could be readily incorporated, or even, given the modest and isotropic change in the elastic properties with (sub)stoichiometry, might be safely ignored.

The examination of the structural phase stability with this stoichiometry yields insight into the nature of the bonding that determines the structure. Consider again the group IVB metals: all crystallize in the tetragonal phase for the dihydride, but all adopt the cubic phase over a range of substoichiometric hydrogen. The single-phase cubic structure, the δ-phase, is present for a $\text{M}$:H ratio of 1.5 to 1.8 for Ti [46], ranging from 1.5 to 1.7 for Zr [47], and appears below  ∼1.7 for Hf [48]. Above this δ-phase, the lattice converts to the ε-phase (fct structure), below this range it becomes a mixed α(hexagonal)-δ phase, essentially hydrogen dissolved in the pure-metal lattice. This trend of conversion to the cubic phases is consistent with a rigid band picture [19], where the band structure of these different hydrides are very similar to one another, and the principle difference is the number of valence electrons that fill that band structure. The $\text{M}{{\text{H}}_{1}}$ (or $\text{MH}_{2}^{0.50}$) would be isoelectronic with the RE cubic dihydrides (the RE atoms having three valence electrons, the IVB atoms four). Within an average-atom picture, the effective substoichiometric HfH$_{2}^{0.75}$ already favors the cubic structure, even before the HfH$_{2}^{0.5}$ that is isovalent. As one adds electrons to the ErH₂ lattice, transmuting to ErH$_{2}^{1.5}$, the crystal distorts to the tetragonal form of the isovalent HfH₂.

Figure 4 plots the computed shear elastic constants, C₄₄ and C_s, as one varies H^z, to fill the valence band structure with more (H$\to {{\text{H}}^{z=1.25}},{{\text{H}}^{z=1.50}}$) and fewer (H$\to {{\text{H}}^{z=0.75}},{{\text{H}}^{z=0.50}}$) electrons. The elastic constants for cubic LaH$_{2}^{z}$, and ErH$_{2}^{z}$ are shown, along with the HfH$_{2}^{z}$ versus the number of electrons filling the valence band structure. Plotting against valence electron count, all three show similar behavior, with negative shear moduli at six valence electrons f.u.⁻¹ signaling the instability of the cubic phase to a distortion. This supports a rigid band picture, where the filling of the band structure rather than the chemical nature of the metal determines the structure. Whether one attributes this to a Jahn–Teller-like distortion due to a high density of states at the Fermi level for the group IVB dibydrides, or disputes that interpretation otherwise through a detailed analysis of the HfH₂ band structure [22], it is unambiguous from comparing these disparate dihydrides that it is the filling within the band structure with six valence electrons that triggers the instability.

Zoom In Zoom Out Reset image size

The parameters crucial for modeling bubble growth and fracture in aging tritides like ErT₂, such as surface energies and mechanical properties [11], are exceedingly difficult to obtain with accuracy from experimental measurements, indicating a first-principles approach as a means to obtain the materials properties [4] is needed to model the bubble evolution, growth and ultimate fracture of the material. The first-principles calculations in this paper provide detailed analysis of the mechanical properties, such as the elastic constants, paying careful attention to the numerical uncertainties that attend the DFT calculations.

The thin film materials are polycrystalline, so the quantities of interest for microstructural modeling of the materials are averaged quantities, such as the bulk shear modulus μ, and the Young's modulus E. For single crystals, these quantities can be obtained straightforwardly from the computed elastic constants [33, 39], but assessment of these quantities in a polycrystalline material adds further numerical uncertainties into the analysis. Analyses by Reuss and Voigt for the shear modulus, ${{G}_{\text{R}}}$ and ${{G}_{\text{V}}}$, place bounds on the bulk shear modulus in a polycrystalline material. Both of these bounds incorporate a numerical uncertainty roughly equivalent to the numerical uncertainty in computing the crystalline shear elastic constants. Conventionally, the 'predicted' bulk shear modulus ${{G}_{\text{H}}}$ is defined as the average of these [33], but in many cases this convention injects additional numerical uncertainties into the 'predicted' shear modulus, as the difference in the ${{G}_{\text{R}}}$ and ${{G}_{\text{V}}}$ bounds can be significant, larger than 10 GPa, for the hydrides.

This shear modulus G, with the bulk modulus B, is used to estimate Young's modulus E:

$$\begin{eqnarray}&&E=\frac{9B{{G}_{\text{H}}}}{3B+{{G}_{\text{H}}}}\end{eqnarray} \tag{ 5 }$$

which, being proportional to ${{G}_{\text{H}}}$, carries roughly the same proportion of numerical uncertainty as G. The consequence is that the numerical calculations for these essential quantities, even obtained from precise DFT calculations, embody numerical uncertainties as large as the rather substantial uncertainties derived from the experimental analysis of these quantities, even discounting the physical uncertainties associated with the accuracy of the DFT approximation.

Table 6 presents the results of the DFT analysis for these bulk mechanical properties, along with their experimental counterparts. The Debye temperature can also be used as an experimental comparison as it can be computed directly from shear moduli [60]. The values in parenthesis are the numerical uncertainties in the DFT predictions. These incorporate uncertainties inherent within the Reuss ${{G}_{\text{R}}}$ and Voigt ${{G}_{\text{H}}}$ bounds from the numerical k-limit uncertainties of the underlying shear elastic constants from which these are computed, and then adds an uncertainty given estimating the shear modulus ${{G}_{\text{H}}}$ as an average midway between these bounds. It is this total uncertainty in ${{G}_{\text{H}}}$ that can be compared meaningfully to the appropriate experimental quantity. Even if the DFT approximation were physically perfectly accurate, these numerical uncertainties from the averaging and k-sampling would remain. The shear modulus across these materials can only be predicted to within 10–20%, roughly half of this uncertainty due to the k-convergence in the computation of the elastic constants, and the other half in the estimation of ${{G}_{\text{H}}}$ within the bounds given by ${{G}_{\text{R}}}$ and ${{G}_{\text{V}}}$.

Table 6. Experimental and computed moduli (in GPa) and Debye temperatures (K).

		Young's modulus	${{G}_{\text{R}}}$—reuss shear bound	${{G}_{\text{V}}}$—voigt shear bound	${{G}_{\text{H}}}$—shear Modulus	Bulk modulus	Debyea temperature (K)
ErH2	Exp.	$148\pm 20$			$60\pm 10$	$97\pm 4$	381
	GGA-PBE	144 (17)	55.2	62.1	59 (7)	89 (1)	389.8 (6.5)
YH2	Exp.	$135\pm 20$			$55\pm 10$	$90\pm 7$	537
	GGA-PBE	136 (15)	52.5	58.6	56 (6)	83 (1)	521 (8)
LaH2	Exp.	$36\pm 6$			$14\pm 3$	$24\pm 3$	243
	GGA-PBE	95 (10)	34.6	41.2	38 (4)	89 (1)	365.5 (6.5)
TiH2	Exp.	$100\pm 15$			$40\pm 7$	$66\pm 5$
(fct)	GGA-PBE	68 (14)	20.8	27.1	24 (5)	143 (1)
ZrH2	Exp.	$175\pm 20$			$70\pm 10$	$115\pm 7$
(fct)	GGA-PBE	113 (24)	35.4	47.7	42 (9)	131 (1)

^aExperimental Debye temperature from [66].

Young's modulus E is, to lowest order, proportional to ${{G}_{\text{H}}}$, see equation (5), and thereby inherits the same proportion of uncertainty as ${{G}_{\text{H}}}$, which amounts to 10–24 GPa, uncertainties comparable to those seen in the experimental measurements. Once more, in the face of such numerical uncertainties, it is pointless to quote ${{G}_{\text{H}}}$ or E to finer than 1 GPa.

The experimental and first-principles-based results for ErH₂ and YH₂ exhibit very good agreement, the difference in ${{G}_{\text{H}}}$ being  ∼1 GPa. We reiterate that this close agreement, while gratifying, is not very meaningful; both the experimentally derived and first-priniciples computed ${{G}_{\text{H}}}$ carry 10 GPa uncertainties. The bulk modulus comparison, where the PBE results are  ∼10% smaller than the experimentally inferred value, is more typical of the kind of agreement that can typically be expected between DFT and experiment. These results both confirm (validate) the accuracy of the DFT for the elastic properties of the dihydrides, and lend greater confidence to the computed values of these mechanical properties, greater confidence than in the rather involved analysis to obtain these values from experiment. The 10% differences in B are expected, and the computations offer refined guidance for microstructural models.

The comparisons between experimental and computational values for LaH₂, TiH₂ and ZrH₂ are much worse. The results for ZrH₂ are almost within agreement, strictly speaking, invoking the full combined uncertainties in both the experimentally-derived values and the DFT-derived values. However, the results for LaH₂ and TiH₂ are clearly well outside of the range of the uncertainties of the experimental analysis and the DFT analysis. These comparisons illustrate the value of using the experimental and computational approaches in concert for such a study. The physical uncertainties in DFT, particularly for B are known to be modest. The PBE is typically prone to be too soft, by 10–20%; this accuracy of DFT for the bulk modulus is confirmed for ErH₂ and YH₂. The likely explanation of the difference between the modeled cubic materials and the experimental results is the difficulty in manufacturing samples of controlled known phase. Both LaH₂ and TiH₂ have complex phase diagrams, and multiple phases in close proximity. Slight variations in manufacturing procedure will produce a variety of phases, possibly even mixed phases, corrupting the analysis of elastic properties of the ideal dihydride material. The problem in ZrH₂ is similar but less extreme: most likely the cubic phase was also in the manufactured material. No x-ray diffraction results were taken, the specific phases or their proportion present in the LaH₂ and TiH₂ samples are not known. What this comparison exposes is the tremendous challenges in obtaining accurate assessment of the mechanical properties from nanoindentation experiments for dihydride films. For these materials, the DFT-derived results can be used with greater confidence in microstructural models than analyses descended from the nanoindentation experiments.

The Debye temperature estimated for ErH₂ and for YH₂ is in extraordinary good agreement with the measured experimental values, as is the measured ${{\Theta}_{\text{D}}}$ for LuH₂ (361 K) [66] to the computed value in table 2. This would appear to validate the relationship between shear moduli and ${{\Theta}_{\text{D}}}$ proposed by Deligoz, et al [60] The discrepancy is larger for LaH₂, experiment measuring 243 K [66] while the DFT result predicts 366 K. At this point, it is not possible to state whether this is a failing in the model, or a failing in the PBE functional, or whether the experimental analysis is in some manner flawed. It would be useful to revisit the experimental measurements of ${{\Theta}_{\text{D}}}$ in light of the sample preparation difficulties encountered in attempting to measure elastic properties in LaH₂, to resolve this discrepancy.