Molecular dynamics studies of fundamental bulk properties of palladium hydrides for hydrogen storage

Solid-state hydrogen storage materials hold the potential to replace current high-pressure storage solutions and accelerate the transition away from fossil fuels.¹ For this technology to become commercially competitive, however, properties of existing hydrogen storage materials must be improved. Improvement of the materials requires a better understanding of the kinetics of hydrogen uptake (hydrogenation) and release (dehydrogenation) processes. Unfortunately, such kinetic processes are complex and are difficult to measure in experiments. Thus, they generally remain poorly understood despite extensive past experimental studies. One effective means to study complex kinetics of (de)hydrogenation processes at realistic time and length scales is through simulations using certain continuum-based techniques, such as phase field models.^2,3 However, the primary challenges with such continuum models are the lack of experimental input parameters and the inability to independently validate each physical relationship within the continuum framework. Higher level computational methods such as molecular dynamics (MD) can both produce inputs and validate assumptions for continuum models without additional experimental information other than the interatomic potential. Since the purpose is to understand the fundamental rules governing a class of materials, we need only consider a representative material towards developing such an MD-continuum multiscale approach.

Palladium (Pd) can efficiently split hydrogen (H) molecules on its surface and has a high solubility and a high diffusivity of dissociated hydrogen atoms in its bulk.⁴ Consequently, hydrogen atoms can easily permeate into Pd to form an interstitial hydride, PdH_x. Thus, Pd is widely used for many hydrogen-related applications, including hydrogen or hydrogen isotope storage, separation membranes, and chemical sensing. Pd is also an important component in hydrogen-based catalytic reactions (e.g., oil reforming, hydrocracking, and hydro-processing).⁵ Among solid-state hydrogen storage materials, palladium hydrides fit within the category of interstitial metal hydrides. Furthermore, the diffusional phase transformation of PdH_x can be considered a fundamental basis for understanding more complicated phase transformation mechanisms in other metal hydrides.⁶ As a result, PdH_x has been the subject of extensive studies, including highly advanced experimental characterization techniques designed to isolate individual phase transformation mechanisms.^7–10 For these reasons, PdH_x can be considered an ideal foundational system for modelling solid-state hydrogen storage reactions.

Within the miscibility gap range, the phase-separating tendency of Pd-H systems creates an equilibrium between α (low hydrogen concentration) and β (high hydrogen concentration) phases. Hydrogenation and dehydrogenation processes proceed through changes of volume fractions of the α and β phases. Once the dissociation/recombination rates between molecular H₂ in the source phase (gas or liquid) and atomic H on the Pd surface are known, (de)hydrogenation processes can be fully determined by specific bulk and non-bulk properties. Bulk properties of interest here include lattice constants, Gibbs free energies of bulk formation, elastic constants, and bulk diffusivities; corresponding non-bulk properties include surface energies, interfacial energies, surface diffusivities, and surface segregation ratios. Together, these properties determine energetic driving forces and kinetic parameters. For example, the coherency strain energy, which is determined by the mismatches in lattice constants and elastic constants between coexisting α and β phases, strongly affects the driving force and, hence, the kinetics of the α ↔ β phase transformations in PdH_x. Similarly, diffusivities directly impact the kinetics, and surface-induced H segregation may modify H surface diffusivities as well as the H₂ ↔ 2H chemical reactivity at the surface. Recently, we have performed molecular dynamics simulations to quantify hydrogen diffusion in bulk PdH_x in the entire reaction space of temperature and composition.¹¹ 

Based on the above discussions, this paper uses PdH_x as a prototypical metal hydride material to study the remaining bulk properties (the non-bulk properties will be addressed separately) with four major objectives: (i) illustrate robust MD methods capable of calculating lattice constants, Gibbs free energies, and elastic constants at any composition and temperature; (ii) perform MD simulations to construct analytical expressions of these bulk properties as functions of temperature and composition for convenience of developing large-scale models (e.g., phase field modelling) to directly simulate (de)hydrogenation; (iii) use MD data to elucidate the experimental observations; and (iv) discuss new physics revealed in the simulations that can help understand solid-state processes in hydrogen storage materials. It is important to note that MD results depend on the interatomic potential which always involve errors. However, our purpose is to develop an understanding of fundamental rules governing a class of materials, which are less dependent on the interatomic potential. In addition, MD can guide development of large-scale models and validate them because when using the MD results as inputs, the large-scale models should produce the same results as MD regardless of the potential. Finally, MD allows the calculations of properties as a continuous function of composition (which cannot be measured from experiments in the two-phase region), and the calculations of statistically averaged properties of large systems in the entire composition and temperature space (which cannot be obtained from quantum mechanical calculations).

To simulate (de)hydrogenation processes in their entirety, any input properties used in phase field or other continuum models must be expressed as a function of both temperature and composition. This imposes some challenges for MD simulations because unlike molecular statics (MS) calculations of these properties at 0 K, finite-temperature MD simulations necessarily involve thermal noise. Time-averaged MD simulations¹² will be explored to resolve this problem. As one example, elastic constants are calculated, using the finite difference equation, as C_ij = δσ_i/δε_j, where δσ_i is a small stress caused by a small applied strain δε_j. The challenge to be overcome here is that the meaningful calculations can only be possible when the thermal noise of the MD-calculated δσ_i is significantly smaller than the δσ_I itself that is already very small. To reduce the thermal noise, we replace snapshot value of δσ_i with a time-averaged δσ_i, which allows us to explore the convergence of C_ij = δσ_i/δε_j with respect to the averaging time.

All simulations were performed using the MD code LAMMPS^13,14 based on the literature Pd-H embedded-atom method (EAM) potential.¹⁵ Due to the small mass of hydrogen atoms (∼1 amu), an extremely small time-step of dt = 0.5 fs was employed to guarantee stable MD simulations except in Sec. V B and  Appendix where effects of larger (artificial) hydrogen mass and time-step are explored.

To cover the entire hydrogen composition range, 11 PdH_x samples with compositions x = 0.0, 0.1, 0.2,…, 1.0 were studied in Secs. III–V to be discussed below. Each sample contained n × n × n unit cells of a face-centered-cubic (fcc) palladium lattice aligned in the three ⟨100⟩ coordinate directions (the value of n for specific simulations will be given below), with the corresponding number of hydrogen atoms occupying initially the octahedral interstitial sites. The system dimensions were relaxed using zero-pressure NPT (constant number of atoms, pressure, and temperature) MD simulations. For clarity in the following discussions, we first present the simulation results and then compare simulations and experiments in dedicated subsequent Sec. IV.

To thoroughly quantify the lattice constants as a function of temperature and composition, 132 zero-pressure NPT MD simulations were performed across a matrix of 11 samples (sample dimension 8 × 8 × 8 unit cells) at 12 temperatures (T = 300.0, 325.0, 350.0,…, 550.0 K). After the first 100 ps were discarded to establish an equilibrium state, lattice constants were derived from the time-averaged system dimensions obtained for the next 2.2 ns of simulation time. Upon fitting the MD results with a polynomial, we obtain the following general equation for the lattice constant a(x,T) (Å) as a function of temperature T (K) and composition x:

In Eq. (1), $a300 K$ and α are the lattice constant (Å) and thermal expansion coefficient at 300 K, respectively. The precision of the parameters in Eq. (1) is carefully chosen to ensure that the calculated lattice constant a(x,T) is precise to the fourth decimal place. The fitted lattice constants as a function of composition are plotted alongside the MD data in Fig. 1 for three selected temperatures: T = 300, 400, and 550 K. It can be seen from Fig. 1 that all MD data points fall on smooth trends with negligible error bars, validating that the time-averaged MD simulations produce well-converged finite temperature lattice constants. In addition, Fig. 1 demonstrates that our fitted Eq. (1) represents the extensive MD data well.

The results shown in Fig. 1 are consistent with our recent discoveries that long time-averaged MD surprisingly produce much more converged results than MS simulations for large and complex systems.^16,17 This can be explained by the fact that unlike MS simulations, which are prone to trapping in metastable atomic configurations, MD simulations directly sample the available configurational space in an unbiased way and can therefore avoid traps with sufficient simulation time. Statistical errors (e.g., thermal noise and numerical truncation) in MD simulations also tend to satisfy normal distributions, which lend themselves well to averaging via increased simulation time. As an example, our previous work¹² indicated that when dislocations were introduced in a big system (∼1 million atoms), the total system energy calculated by two independent MS configurations for a given composition differed by ∼100 eV (∼0.0001 eV/atom). On the other hand, thermal averaging from different initial configurations gave time-averaged MD results that were within 1 eV. For the simple property of a small system explored here, MS is expected to also produce highly converged results. To verify this, MS simulations were performed to calculate the lattice constants at various compositions, and the results are included in Fig. 1. Clearly, the MS data also fall on a smooth trend. Note that the MS simulations necessarily use randomly initialized spatial population for hydrogen atoms.

One motivation to use relatively small systems (8 × 8 × 8 unit cells) is to prevent phase separation so that we can access continuous composition and temperature space that may be not under equilibrium and therefore inaccessible in experiments. This is because in small systems, any phase separation will cause large interface boundary areas, and defects that can relieve misfit strains between phases (e.g., misfit dislocations) are less likely form. To verify that our system did not involve significant hydrogen segregation, we explored the final configuration obtained for the PdH_0.4 sample undergoing a simulation at 600 K. The y-z cross section of the sample was divided into 2 × 2 columns, and hydrogen density along each column was calculated as a function of x. The results are shown in Fig. 2. Figure 2 confirms that the hydrogen distribution is relatively uniform. The uniform hydrogen population is further supported by our recent investigation into the dynamics of H in small system sizes of PdH_x,¹¹ which demonstrated a diffusive behavior that could be well described by a single-barrier random-walk process for compositions x < 0.7.

Next, we calculate the Gibbs free energy of formation of bulk PdH_x. First, we provide an analytical definition of our Gibbs free energy of formation. Considering that each PdH_x chemical unit contains 1 + x atoms, the chemical reaction to form “one atom” of PdH_x from the two reference states Pd and PdH can be written as

The heat of formation ΔH_f in the per-atom unit is defined as the enthalpy change associated with the reaction in Eq. (2)

where E_PdHx, E_Pd, and E_PdH are the energies of PdH_x, Pd, and PdH, respectively, and X is the mole fraction of hydrogen. Here, all energies use the per-atom unit. For instance, the total energy for the chemical formula PdH equals 2E_PdH. Note that enthalpies H can be represented by energies E because our simulations were conducted at zero pressure. Finally, mole fraction X and H/Pd ratio x are interchangeable through relations $X=x1+x$ and $x=X1−X$⁠.

Using the same procedures as described above for determining lattice constants, time-averaged $EPdHx$⁠, $EPd$⁠, and $EPdH$ values were calculated over the 2.2 ns period in the same 132 MD simulations. The results were then used to calculate the heat of formation using Eq. (3). By fitting the MD data of heat of formation to the Redlich-Kister polynomial^18–20 with the polynomial coefficients expressed as parabolic functions of temperature, we obtain a general equation of heat of formation ΔH_f(X,T) as a function of temperature and composition as (in the unit of eV/atom)

Again, the precision of the parameters is carefully chosen so that the calculated ΔH_f(x,T) is precise to the fourth decimal place. The polynomial expression was chosen for convenience of fitting, with the order chosen to ensure adequate reproduction of the MD results without overfitting. The fitted curves and the MD data of ΔH_f(X,T) as a function of composition are compared in Fig. 3(a) at two selected temperatures (300 and 400 K). Like the lattice constant case, we also include MS data in Fig. 3(a). Clearly, the MD ΔH_f(X,T) data carry a negligible statistical sampling error with smoothness comparable to the MS results, and our fitted Eq. (4) reliably captures the simulated MD data at different compositions and temperatures.

Due to the short-range ordering, accurate calculations of the configurational entropy of PdH_x requires Monte Carlo simulations.²¹ However, our systems are rather uniform as can be seen from Fig. 2. Hence, we calculated an approximate configurational entropy based on the assumption that hydrogen atoms (or equivalently vacancies) in the fcc octahedral sublattice are randomly distributed. Under this assumption, the entropy of formation in a per-atom unit can be written as¹⁵ 

where k_B is the Boltzmann constant. The Gibbs free energy of formation is then

The Gibbs free energy of formation curves calculated from Eqs. (4) to (6) are shown in Fig. 3(b) for three selected temperatures: 300, 500, and 700 K. Figure 3(b) confirms that the Pd-H potential employed in this work produces the double-well type Gibbs free energy, resulting in phase separation known for the system.²² Specifically, a common tangent intersects with the 300 K curve at approximately two compositions X_α = 0.03 and X_β = 0.5. The comparison of these predicted equilibrium compositions with corresponding experimental values will be further discussed below.

Finally, we point that Eq. (5) ignores the vibrational entropy contribution. To justify this, we explicitly computed the vibrational entropy of 11 samples (x = 0.0, 0.1,…, 1.0) at 300 K, with details given in the supplementary material. We found that the vibrational entropy of Pd is well described by a linear function of composition (almost a constant). This indicates that the vibrational entropy of formation, which is defined as the difference in vibrational entropies between hydride PdH_x and a linear combination of Pd and PdH, is negligible.

Conventionally, direct calculations of elastic constants using C_ij = ∂σ_i/∂ε_j are performed at zero temperature based on MS simulations. Calculations of finite-temperature elastic constants require MD simulations. Normally, such MD calculations are achieved through fluctuation theory^23,24 or energy vs. strain curves^16,25 rather than from explicit calculation of C_ij = ∂σ_i/∂ε_j. Direct calculations are sometimes desired, since they not only lead to straightforward determination of all elastic constants C_ij(i, j = 1, 2,…, 6) but also reveal insights that are difficult to observe using alternative approaches. Such calculations, however, may seem formidable at first sight because the thermal noise associated with the MD simulations could make it difficult to accurately calculate ∂σ_i/∂ε_j using the finite difference method. Encouraged by our recent success in using time-averaged MD simulations to calculate finite-temperature Young's moduli,¹⁶ here we explore the feasibility to directly calculate highly converged finite temperature elastic constants C_ij (i, j = 1, 2,…, 6) for PdH_x. Two direct methods are studied. The first one uses a steady loading normally applied in mechanical testing experiments, and the second one uses a vibrational loading normally applied in the more common ultrasonic testing experiments.

The key to calculate elastic constants from C_ij = δσ_i/δε_j is to reliably extract the small stress response δσ_i at a small imposed strain δε_j. Within MD simulations, sufficient time averaging is required to reduce the thermal noise of δσ_i to a negligible level. To maintain the imposed strain, the MD simulations should be performed under an NVT condition (constant number of atoms, volume, and temperature). To ensure the equilibrium dimension of the crystal at finite temperature, a trial MD simulation using a zero-pressure NPT ensemble was first performed for 100 ns. After the first 20 ns was discarded, the time-averaged lattice volume was calculated for the final 80 ns. A new crystal was then created based on the time-averaged dimensions, thereby implicitly accounting for thermal expansion. After further relaxation for 100 ns in an NVT MD simulation to ensure equilibration, this carefully prepared equilibrium crystal was then used to calculate finite-temperature elastic constants.

Starting from the equilibrium crystal, positive and negative small strains of the j-th component, ±ε_j = 10⁻⁴, were separately applied to the system. An MD simulation using an NVT ensemble was performed for 100 ns to relax both the positively and negatively strained systems. After discarding the first 20 ns, time-averaged reduced stresses ⟨σ_i⟩ (i = 1, 2,…, 6) were calculated for the remaining 80 ns. The MD elastic constants C_ij,MD (i = 1, 2,…, 6) were then calculated using a finite-difference scheme

By repeating the same process for all j = 1, 2,…, 6, we determine all elastic constants. According to the symmetry of the cubic crystal, only C_ij (i, j = 1–3) and C_ii (i = 4–6) are non-zero tensor components.²⁶ Further, these non-zero components satisfy C₁₁ = C₂₂ = C₃₃, C₄₄ = C₅₅ = C₆₆, and C₁₂ = C₁₃ = C₂₁ = C₂₃ = C₃₁ = C₃₂. However, the elastic constants C_ij,MD obtained from MD simulations do not strictly satisfy these relations because of their statistical nature. Hence, we recast the elastic constants using the following equations:

The recast C_ij enforce the cubic relations described above, where only C₁₁, C₁₂, and C₄₄ are independent. These three independent elastic constants can be converted to alternative independent elastic constants: Bulk modulus B = (C₁₁ + 2C₁₂)/3, shear modulus C′ = (C₁₁ − C₁₂)/2, and shear modulus C₄₄. To provide objective measures on the quality of the elastic constant calculations, we define four standard deviation parameters ξ_i (i = 1,…, 4) of the raw MD elastic constants from the four cubic constraints

The framework described above was applied to calculate elastic constants for 11 sample compositions x = 0.0, 0.1, 0.2,…, 1.0 (sample dimension 4 × 4 × 4 unit cells), each at seven temperatures (10, 100.0, 200.0,…, 600.0 K). Here, smaller system sizes help reduce the computational cost. They also further guarantee a single uniform phase without phase separation. The calculated results of bulk modulus B and shear moduli C′ and C₄₄ are shown in Figs. 4(a)–4(c), respectively, using black symbols, and the four error parameters ξ_i (i = 1,…, 4) for all the samples displayed in Figs. 4(a)–4(c) are shown in Fig. 4(d). For comparison, the same elastic constants obtained from 0 K MS simulations are included in the figures using black stars. The MD elastic constants C (C = B, C′ and C₄₄) are fitted to the following equation:

The fitted parameters c₁, α₁, c₂, α₂, β₂, c₃, α₃, γ₃, β₃, and x₀ are shown in Table I, and the fitted curves are shown as black lines in Figs. 4(a)–4(c). We also include literature experimental data^27–29 using colored symbols, although we defer discussion of the comparison to the experimental values to dedicated Sec. IV below.

Multiple important conclusions can be drawn from Fig. 4. First, it can be seen from Fig. 4(d) that in general, the errors of our calculations based on Eqs. (9)–(12) are negligibly small. This strongly validates that time-averaged MD simulations can indeed provide a robust method to calculate finite temperature elastic constants. Second, certain data points at lower temperatures, specifically those at 100 K and 200 K, are excluded in Figs. 4(a)–4(c). The reason that those data points are omitted is because they are associated with much larger errors [The errors are beyond the range of Fig. 4(d) if included.] than the results at the other temperatures. To confirm that these errors are not due to inadequate statistical sampling, we repeated the same extensive simulations in  Appendix using different random number seeds.  Appendix also explores the use of artificially increased hydrogen mass to accelerate simulations. Almost the same results are seen in  Appendix, where again the cubic deviation error parameters remain excessive only at 100 K and 200 K. Thus, these errors arise from physical reasons that we do not yet understand completely. As plausible explanations, although no phase separation is observed in our simulations due to the small computational system size, PdH_x at low temperatures has a large thermodynamic tendency towards phase separation, as can be seen from Fig. 3(b). This tendency may cause the system to violate the cubic elastic constant relations. The tendency for phase separation may be kinetically inhibited when the temperature is further reduced. Therefore, it is not surprising that it is difficult to obtain cubic elastic constants at 100 K and 200 K but not at lower or higher temperatures. Third, it can be seen from Figs. 4(a) to 4(c) that while MD data points seem to exhibit minima whereas Eq. (13) monotonically decreases with composition, the deviation of Eq. (13) from the MD data is not significant. This means that the use of Eq. (13) in continuum models will not lead to significant errors as compared to MD simulations. Finally, Figs. 4(a)–4(c) also reveal a significant and intriguing temperature effect: The PdH_x elastic constants monotonically decrease with composition at temperatures T ≥ 100 K, whereas they are not monotonic with composition at 0 and 10 K. A related phenomenon is that elastic constants at T ≥ 100 K are much smaller than those at 0 K for high hydrogen compositions.

Vibrational loading is important to explore because most elastic constant measurements are performed using ultrasonic methods rather than direct mechanical testing.^27–29 Under a sinusoidal vibrational loading condition, the stress σ_i caused by a sinusoidal strain ε_j (all other strain components are assumed to be zero) can be expressed by Hook's law as

where $σi+ι$ is an instantaneous stress measured from MD that necessarily involves thermal noise ι; ε_j,0 and Δε_j are the average (residual) and the amplitude of the imposed strain, respectively, and τ is the period of the strain oscillation. To prevent the unnecessary cancellation of positive and negative stresses/strains, Eq. (14) is squared before taking the time average. Under the condition that the total simulation time t is a multiplication of τ, we can derive the following equation:

Equation (15) means that elastic constants can be evaluated either from the slope of the linear $〈(σi+ι)2〉$ vs. $Δεj2$ relationship at a given residual strain ε_j,0 or from the slope of the linear $〈(σi+ι)2〉$ vs. $εj,02$ relationship at a given amplitude $Δεj2$⁠. This is important because the calculations of these slopes do not require precise information about thermal noise ι and residual strain ε_j,0, Note that $〈(σi+ι)2〉$ converges to $〈σi2〉+〈2σi⋅ι〉+〈ι2〉=〈σi2〉+〈ι2〉$ when t → ∞.

The same equilibrium crystals used for the steady loading were used for dynamical loading. To ensure converged results for the more challenging dynamical loading, total simulation time was increased to 800 ns. After the first 80 ns was discarded, $〈(σi+ι)2〉$ was calculated for the remaining 720 ns. To keep the computation tractable, we artificially increased the hydrogen mass in the vibrational loading simulations to 50 amu. This allowed use of a much larger time step (dt = 4 fs rather than 0.5 fs), meaning that we can increase the total simulated time by almost an order of magnitude with the same computing resources. We demonstrate in  Appendix that the results are not sensitive to this artificial increase in mass. To further verify that our dynamical loading simulations with increased mass and time step produce consistent results as shown in Fig. 4, seven test simulations were performed at seven residual strains ε_4,0 = −2.0 × 10⁻⁴, −1.5 × 10⁻⁴, −1.0 × 10⁻⁴, 0.0, 1.0 × 10⁻⁴, 1.5 × 10⁻⁴, 2.0 × 10⁻⁴, a fixed amplitude of Δε₄ = 0, a fixed composition of x = 0.0, and a fixed temperature of 500 K. Because Δε₄ = 0, we are essentially emulating steady loading, which can also be viewed as starting with a strain ε_4,0 and then applying the sinusoidal loading with the period τ → ∞. The resulting ⟨(σ₄ + ι)²⟩ vs. (ε_4,0)² data are shown in Fig. 5. An ideally linear ⟨(σ₄ + ι)²⟩ vs. (ε_4,0)² plot is obtained. Hence, we can reliably determine C₄₄ from the slope, which gives C₄₄ = 39.9 GPa. This is in a very good agreement with a C₄₄ of 41.8 GPa determined in Fig. 4(c) under the steady loading condition, thereby directly validating our vibrational loading methods.

Having justified the method, MD simulations are performed to calculate the vibrational loading elastic constant C₄₄ at two compositions of x = 0.0 and x = 0.7 and one temperature of T = 500 K. Because these simulations are much more expensive, the scope of the calculations was reduced so that we could focus on the physical significance of the dynamical loading.

For each composition, 56 simulations consisting of a matrix of eight periods τ = 72, 80, 90, 120, 144, 180, 240, and 360 ns and seven strain vibrational amplitudes of Δε₄ = 2.0 × 10⁻⁴, 2.5 × 10⁻⁴, 3.0 × 10⁻⁴, 3.5 × 10⁻⁴, 4.0× 10⁻⁴, 4.5 × 10⁻⁴, and 5.0 × 10⁻⁴ were performed for 800 ns. Although the nominal value of the residual strain ε_4,0 was set to zero, ε_4,0 is unlikely to be exactly zero during NVT simulations. Thus, it was treated as a fitting parameter. After the first 80 ns was discarded, $〈(σi+ι)2〉$ was calculated for the remaining 720 ns (which is a multiple of the period τ in all cases). As an example, the $〈(σi+ι)2〉$ vs. (Δε₄)² data obtained at two oscillation periods of τ = 72 and 360 ns are shown in Fig. 6(a) for the case of pure Pd (i.e., x = 0). Linear relationships are obtained between $〈(σi+ι)2〉$ and (Δε₄)², confirming the validity of Eq. (15) for calculating elastic constants under vibrational loading modes. Based on the slopes of the $〈(σi+ι)2〉$ vs. (Δε₄)² relationships, the elastic constant C₄₄ was calculated for different vibrational periods τ, and the results are shown in Fig. 6(b). The computed C₄₄ for pure Pd varies within a range between 40 and 75 GPa for the oscillation periods we considered. Because the relationships shown in Fig. 6(a) are ideally linear, we conclude that the variation of C₄₄ with τ is a physical response rather than statistical error. Intriguingly, significant variation occurs in the shorter vibration period regime as one can see from Fig. 6(b). Therefore, this variation is likely to be caused by the mismatch between the loading oscillation period and the mechanical response time.

Similar ⟨(σ₄ + ι)²⟩ vs. (Δε₄)² data are obtained for the higher composition of x = 0.7, as shown in Fig. 7. In contrast to the case of pure Pd, these data do not exhibit a linear relationship. This means that constant elastic constants cannot be deduced from this analysis, at least on the time scales we employed. Again, because our long time-averaged MD simulations should not cause significant statistical errors, the results shown in Fig. 7 corroborate Fig. 4 that the presence of hydrogen mechanistically changes the elastic response at finite temperatures. We note that significant H diffusion occurs under our simulated conditions; hence, the non-linear behavior seen in Fig. 7 may be caused by the modulated stress response due to hydrogen diffusion. We will return to this important point in the discussions below, including possible implications for interpreting the experimental data in Fig. 4.

The 0 K lattice constant predicted by the interatomic potential is known to produce accurately the corresponding experimental value for Pd.^15,30 Our own Rietveld refinements for the XRD pattern measured on a bulk Pd powder indicate a lattice constant of 3.895(2) Å at 300 K, which is also consistent with the calculated value shown in Fig. 1. In addition, the room-temperature palladium thermal expansion coefficient determined in Eq. (1), α ≈ 14 × 10⁻⁶K⁻¹, is quite close to the value of α ≈ 12 × 10⁻⁶K⁻¹ determined in experiments.³¹ On the other hand, Fig. 1 indicates a PdH lattice constant of above 4.3 Å at 300 K, which is significantly larger than the experimental value of 4.09 Å at 77 K.³² The actual deviation may not be that big because the equilibrium composition of β phase achieved in experiments is likely to be lower than PdH [X_β = 0.38 or x_β = 0.61 (Ref. 22)].

The experimental equilibrium compositions of α and β phases at 300 K are approximately X_α = 0.02 and X_β = 0.38,²² as compared to our predicted values of X_α = 0.03 and X_β = 0.5 seen in Fig. 3. Consequently, the potential overestimates the composition of the β phase.^10,33 Note that with model errors on PdH lattice constant and equilibrium β composition established, we can more clearly understand the general trends and derive general rules governing a class of materials, namely interstitial hydrides.

Figure 4 indicates that our calculated 0 K bulk modulus B and shear modulus C^′ match well with the literature ultrasonic measurements.^27–29 C₄₄ deviates more from the experimental data^27–29 in terms of value, but its trends with composition are similar to B and C^′. These elastic constants initially decrease and then increase with composition until they reach a maximum at a high composition near x ∼ 0.6.^27–29 In contrast, the elastic constants calculated at temperatures T ≥ 300 K nearly monotonically decrease with composition. As a result, the elastic constants calculated at these temperatures are significantly lower than the ultrasonic experimental values when hydrogen composition is high even though elastic constants of Pd (B and C^′) still match well with experiments. This is extremely interesting because this discrepancy cannot be attributed to inaccuracies in the interatomic potential given the good agreement at 0 K (especially because interatomic potentials are usually fully defined by the 0 K energy profiles). Although the overestimated PdH lattice constants by the potential may result in reduced elastic constants, this explanation can be excluded because it should also impact the 0 K data. Note that previous MD simulations using a completely different interatomic potential indicates a similar decrease in bulk modulus with the hydrogen concentration.³⁴ We would also reemphasize that by sufficiently sampling the associated phases, our MD simulations do not create additional errors beyond the interatomic potential. Evidences verifying this include the near zero statistical error shown in Figs. 1 and 3(a), the near zero deviation of elastic constants from the physical cubic relations as shown in Fig. 4(d) (note that if anything goes wrong, this error would go out of chart), the perfect match between data points and function lines for the Arrhenius plots,¹¹ and an excellent agreement between two independent sets of calculations of elastic constants as shown in Figs. 4 and  Appendix. The cubic elastic constants also verify that our relatively small system dimensions successfully prevented phase separation. Note that the MD and experimental discrepancies of finite-temperature elastic constants are found not just in the magnitudes but also in their trends. Different trends suggest different underlying physics.

The experimental data shown in Fig. 4 are all from ultrasonic measurements.^27–29 The calculated values shown in Fig. 4, however, pertain to mechanical testing with a steady loading. We therefore perform a careful analysis on the literature Pd and PdH_0.6 mechanical testing experimental data.^35,36 We find that the experimental Young's modulus of PdH_0.6 measured from mechanical testing is likely to be significantly lower than that from ultrasonic experiments, although the coarse scale used in the mechanical testing^34,35 may introduce a systematic error. Given that the stresses within the materials under operation are more appropriately described by mechanical testing than by ultrasonic conditions, our finding of a discrepancy between ultrasonic experiments and simulations is important and warrants further careful mechanical testing measurements on elastic properties of palladium hydrides.

The finding that for some special materials, elastic constants measured from ultrasonic experiments can be much higher than those measured from mechanical testing experiments is certainly not new. For instance, extensive past experiments have indicated that for mineral materials, elastic constants measured from vibrational loading can be significantly higher than static loading (even by one order of magnitude).^37,38 PdH_x is a special material. In principle, if the atomic bonds are maintained and only undergo elastic distortions (such as stretching and twisting) in response to mechanical loads, then elastic constants have the traditional meanings. However, PdH_x can be viewed as composed of a conventional solid (the fcc Pd crystal) and a liquid-like H phase trapped in the Pd lattice interstices. In such a system, significant hydrogen migration and repopulation occur. As a result, Pd-H bonds are continuously destroyed and reformed. This hydrogen repopulation can be inhibited either at low temperatures, low composition, or when the required diffusion time cannot be achieved within the short period of a high-frequency vibrational loading. From this observation, it is reasonable that the elastic constants measured from ultrasonic methods may not correspond to the ones measured from mechanical testing for hydrogen-rich PdH_x. The severe interference of hydrogen diffusion within the elastic response is also suggested by the lack of a clear linear relationship in Fig. 7.

In view of the explanations described above, it is worth pointing out that vibrational loading is also affected by the diffusive properties of fully hydrogenated PdH, where significant hydrogen diffusion is also observed within our calculations. This is surprising at first sight because diffusion can only occur through vacancy migration mechanisms. Nonetheless, the apparently fully occupied PdH still contains high concentrations of hydrogen vacancies. In fact, an fcc Pd lattice has two tetrahedral interstitial sublattices in addition to the octahedral interstitial sublattice. Occupying only the octahedral sites will result in a rock salt PdH crystal structure. Alternatively, occupying the first and then the second tetrahedral sublattice of a Pd crystal would give correspondingly a zinc-blende and a CaF₂ structure. A geometrically saturated structure that is free of vacancies can only be obtained when all the three sublattices are fully occupied. However, PdH₃ is thermodynamically unstable and may never occur, meaning that significant numbers of vacant sites remain throughout the hydrogenation process.

Time-averaged MD simulations have been applied to generate high-quality, finite-temperature bulk thermomechanical properties for all compositions of PdH_x. Our study provides not only a more complete understanding of atomistic mechanisms of bulk behaviors during hydrogenation and dehydrogenation processes, but also necessary material parameters for informing mesoscale computational models such as phase-field models. The following achievements are useful both for future simulations and for improving understanding of hydrogen storage applications:

1. We have demonstrated the validity of robust time-averaged MD methods for computing highly converged finite-temperature elastic constants C_ij (i, j = 1,…, 6) under both steady and vibrational loading conditions;

2. We have fitted MD data for lattice constants, Gibbs free energies, and elastic constants to analytical equations as functions of temperature and composition. We have also demonstrated that these fitted functions capture the MD data extremely well. In several instances, significant variations with temperature and composition are observed, suggesting that approximations based on static values cannot accurately capture full hydrogenation/dehydrogenation behavior;

3. Our calculated 0 K elastic constants match reasonably well with the ultrasonic experimental data. However, we found a significant discrepancy between simulations and ultrasonic experiments for elastic constants at temperatures T ≥ 300 K. This discrepancy elucidates well the literature mechanical testing measurements of Young's modulus that appear to differ from the ultrasonic data. Our studies, therefore, indicate that the ultrasonic and mechanical testing of elastic constants of PdH_x may be inconsistent. We hope that this intriguing observation may stimulate further experimental validation studies.

See supplementary material for a complete list of figures showing phonon density of states and the related vibrational entropy as a function of hydrogen composition.

Sandia National Laboratories is a multi-mission laboratory managed and operated by the National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy's National Nuclear Security Administration under Contract No. DE-NA-0003525. Lawrence Livermore National Laboratory is operated for the U.S. Department of Energy under Contract No. DE-AC52-07NA27344. The authors gratefully acknowledge research support from the U.S. Department of Energy, Office of Energy Efficiency and Renewable Energy, Fuel Cell Technologies Office, under Contract Nos. DE-AC04-94AL85000 and DE-AC52-07NA27344. Discussions with J. A. Zimmerman, N. Bartelt, R. Sills, D. Robinson, and D. Cowgill are also highly appreciated.

The views expressed in the article do not necessarily represent the views of the U.S. Department of Energy or the United States Government.

To provide a reference for future researchers who do not have access to supercomputing resources, we repeated the calculations for elastic constants under steady loading, Fig. 4, after artificially increasing the hydrogen atomic mass from 1.01 amu to 50 amu. This allowed us to increase the time step from dt = 0.5 fs to dt = 4 fs. The results of the repeated simulations are shown in Figs. 8(a)–8(c), respectively, for bulk modulus B and shear moduli C′ and C₄₄. Figures 8(a)–8(c) are almost identical to Figs. 4(a)–4(c). For additional certainty, we overlay the 500 K bulk modulus obtained from the two masses in Fig. 8(d). Indeed, Fig. 8(d) confirms that elastic constants are not sensitively affected by the atomic mass. Interestingly, the same set of data points at 100 and 200 K is omitted from Figs. 4 and 8 because they violate the cubic elastic constant relations in both cases. Hence, the reproducibility of this set of data in two independent sets of simulations strongly validates our conclusion that the violation of the cubic relations is not due to statistical reasons but is more likely due to the mechanistic reasons described above.

The lack of difference between Figs. 4 and 8 is also interesting because it contradicts with experimental finding that elastic constants are sensitive to hydrogen isotopes.^27,39 This suggests that the effects of isotopes on elastic constants may be beyond mass and that classical MD methods with NVT ensemble may not capture these effects. Future studies are needed to clarify this problem.