1. Introduction
Plutonium oxide (PuO
) is a stable ceramic that forms spontaneously from the metal in the presence of atmospheric oxygen or water [
1]. PuO
has a number of industrial uses, including as a component in mixed oxide fuels for nuclear power reactors [
2] and as a fuel source for deep-space spacecrafts [
3]. Similar to many transition metals, the oxide layer on plutonium-based materials acts as a protection against corrosion due to atmospheric gases [
4]. In particular, PuO
is thought to help prevent nucleation of hydriding due to the poor solubility of hydrogen in the lattice [
5]. Ab initio molecular dynamics studies have observed that hydrogen diffusion could be slowed by the formation of O–H bonds within the bulk [
6]. Despite this, plutonium metal is still highly sensitive to the presence of adsorbed hydrogen, where the Pu-hydride product can form spontaneously [
7,
8] and strongly exothermically [
9,
10,
11]. Hydriding studies would benefit from atomic-scale calculations to help elucidate the initial stages of this process, including hydrogen absorption and diffusion in the bulk.
In this regard, quantum mechanical calculations using Density Functional Theory (DFT) hold promise as a method to help determine many observables for the H/PuO
system. DFT has been shown to yield accurate results for many bulk PuO
physical and mechanical properties, including the bulk modulus and band gap [
12,
13]. However, many challenges with performing these calculations remain. To the best of our knowledge, there are no systematic studies to determine the best choice of exchange-correlation functional for hydrogen-containing PuO
systems. In addition, use of an onsite Hubbard-like correction (e.g., GGA + U) [
14] combined with calculation of spin effects [
15] is required to yield an accurate lattice constant and band gap. Collinear spin polarization calculations tend to be more computationally efficient [
16], but neglect spin–orbit coupling and relativistic effects that are known to be significant for plutonium [
13]. Both GGA + U and spin–orbit coupling are generally computationally cumbersome due to their creation of numerous local minima in the electron density [
17,
18,
19], which make it difficult to find the ground-state of the system. The magnetic ground-state of PuO
is itself somewhat unknown, with both experiments and DFT calculations yielding a wide range of possible spin states [
12,
13], and the search for the true ground-state could require placing constraints on the orbital occupations [
18,
19]. In contrast, the 1s orbitals from the hydrogen have negligible overlap with the plutonium 5f electrons [
20], indicating that these effects could be less significant for the H/PuO
systems of interest here.
These issues make it difficult to determine basic kinetic properties of the initial hydriding reactions, such as hydrogen diffusion activation barriers and hopping rates. In particular, quantum nuclear vibrational effects (i.e., zero-point energy effects and/or tunneling) could prove large for these systems. Quantum zero-point energies can cause computed hydrogen surface adsorption energies in
-uranium to shift by as much as 30% [
21]. Quantum tunneling has been shown to enhance hydrogen diffusion on transition metal surfaces [
22] and on polycrystalline water ices by several orders of magnitude at low temperatures [
23]. Recent path integral simulations have shown that quantum nuclear vibrational effects can have similar results for proton diffusion in high-pressure transition metal clathrates [
24]. DFT calculations on plutonium containing materials in general are far too computationally intensive to allow for these types of calculations, which can require millions of force evaluations per diffusion pathway. Thus, plutonium hydriding studies could greatly benefit from determination of the appropriate level of DFT theory for computing hydrogen/substrate interactions, as well as a method for placing a bound on the importance of quantum nuclear vibrational effects.
In this work, we attempt to reconcile these issues by first assessing different levels of theory in order to determine whether more computationally efficient DFT calculations could yield satisfactory results for hydrogen absorbed in PuO. We report on results for bulk PuO from several different functionals and spin states, using GGA + U with either spin polarization or spin–orbit coupling in order to assess the effects of each. We then compute hydrogen absorption energies in the dilute limit and determine the relative importance of GGA + U and treatment of spin on hydrogen-PuO interactions. Finally, we use a newly developed semi-empirical quantum Density Functional Tight Binding (DFTB) model to seed DFT calculations of the hydrogen diffusion activation energies. We use these values in a simple quantum double well formalism in order to place a bound on the quantum tunneling effects on single hydrogen diffusion in a PuO lattice. Our results serve as an initial estimate for the importance of considering hydrogen quantum tunneling in these systems, which could play a significant role that requires further consideration.
2. Computational Details
PuO
has been determined experimentally to have a cubic lattice with
symmetry (
Figure 1). The high symmetry of the crystal yields two possible absorption sites for hydrogen: (1) a lower energy “interstitial” site, where an O–H bond is formed and is oriented towards the center of a nearby Pu-Pu axis (labeled
int), and (2) a higher energy “octahedral” site, where the hydrogen ion sits between two Pu atoms (labeled
oct). We focus our efforts specifically on the energetics of these interactions, as well as the the kinetics of proton hops between these sites.
All DFT calculations reported in this work were performed with the Vienna ab initio Simulation Package (VASP) [
25,
26,
27] using projector-augmented wave function (PAW) pseudopotentials [
28,
29]. We have surveyed three different exchange-correlation functionals: the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation (GGA) functional [
30], the reformulation of PBE for solids and surfaces (PBESol) [
31], and the Strongly Constrained and Appropriately Normed meta-GGA semilocal density functional (SCAN) [
32]. Similar to previous studies [
20], we have used the D3 dispersion energy correction [
33] in all calculations discussed here, regardless of treatment of spin or choice of functional.
Lattice optimizations were performed on PuO
alone, using a four formula unit (12 atom) cubic unit cell with a planewave cutoff of 500 eV, Mohkhorst-Pack [
34] k-point mesh of
, energy convergence criteria of
eV and force convergence criteria of
eV/Å, which yielded system pressures far below 1 GPa. Tests with PBE using a planewave cutoff of 600 eV, k-point mesh of
, energy criteria of
eV, and force criteria of
eV/Å indicated convergence of energies within
eV per formula unit and lattice constants within 0.01 Å. All lattice optimizations were performed anisotropically, i.e., without any constraints on the final symmetry of the supercell. In general, this resulted in off-diagonal lattice vector components on the order of 10
Å. Hence, all optimized supercells were assumed to exhibit at least orthorhombic symmetry. Hydrogen absorption calculations were performed for a 32 formula (96 atom) unit simulation cell, using a k-point mesh of
, energy criteria of
eV, and force criteria of
eV/Å.
Use of a linear response approach [
35] to compute
with PBE and PBESol yielded a value of 3.99 eV for both functionals. However, this value resulted in lattice constants and band gaps below experiment for all magnetic states discussed here. Hence, we have used a value of
eV for our calculations with
eV and
eV, based on results from previous work [
13] and independent of choice of functional.
Given the number of different DFT calculations discussed in this work, we adopt a nomenclature to differentiate between different features. First, we consider different collinear and noncollinear periodic arrangements of spin moments, i.e., spin geometries (
Table 1). Diamagnetic states (no magnetic moment) are labeled DM, ferromagnetic states are FM, and antiferromagnetic states are AF. FM states are further labeled (001) for collinear and (111) for three-dimensional noncollinear. For the AF states, we consider one-dimensional collinear geometries which we label ‘1k’, and three-dimensional noncollinear geometries that we label ‘3k’. Each of these can be further differentiated with either longitudinal (long) or transverse (trans) domains, depending on the symmetry of the spin geometry [
13]. For each DFT result, we first state the type of magnetism present, DM, FM, or AF, followed by the spin configuration in parentheses, i.e., (001) or (111) for FM systems, and (1klong), (1ktrans), or (3klong) for AF systems. Thus, the following magnetic states are considered in this work: DM, FM (001), FM (111), AF (1klong), AF (1ktrans), AF (3ktrans).
We have performed calculations using both spin polarization (SP) and spin–orbit coupling (SOC). We note that in the absence of SOC the energy of the system in principle will not depend on the direction of the magnetic moment since the spin moments are not coupled to the orbital angular momenta. Hence, we would expect the total energy to be independent of choice of spin quantization axis. However, the final energy, geometry, and lattice vectors achieved during optimization can still be sensitive to choice of initial spin geometry [
19], where we have observed changes in the magnitude and direction of the spin moments to occur regularly during our calculations. This effect can be pronounced when asymmetries are introduced into the lattice, such as a hydrogen point defect. As a result, we have conducted noncollinear SP calculations as a systematic way of providing additional sampling of SP spin states. The majority of our magnetic calculations have been performed at the SP level. Hence, results from calculations with spin interactions should be assumed to be from SP unless labeled with SOC. Similarly, calculations performed with GGA + U are labeled as such, and the absence of this label implies no GGA + U in that result.
DFTB calculations with self-consistent charges (SCC) were performed using the DFTB+ code [
36], with a many-body repulsive energy (
) determined from the Chebyshev Interaction Model for Efficient Simulation (ChIMES) [
37,
38,
39,
40]. Briefly, the DFTB total energy is derived from an expansion of the Kohn-Sham DFT energy to second order in charge flucuations [
41]. This yields the following total energy expression:
Here,
corresponds to the band structure energy computed from diagonalization of the DFTB Hamiltonian, and
is the self-consistent charge transfer term.
is the empirical “repulsion” term determined from ChIMES, and
is the dispersion energy, determined by the D3 method in our study. Additional details regarding our DFTB parameterization are discussed in the
Supplementary Information.
ChIMES repulsive energies are determined through fitting linear combinations of Chebyshev polynomials to DFT reference data [
38,
39,
40,
42]. We have created an
by fitting a ChIMES potential of up to two-body polynomial order of 12 and three-body polynomial order of 8 order to 75 DFT-MD snapshots for bulk PuO
. These were run at 600 K and at both the experimental density and under 5% compression using DM calculations and no GGA + U interaction in order to create a sufficient data set for our potential. We also have included two PuH
DFT-MD snapshots in order to add some sampling of Pu-H interactions, computed at 1000 K at the experimental density and under 5% compression, and also without spin or GGA + U. These minimal DFT calculations allowed us to sample a large enough part of the phase space with sufficient accuracy for the purposes of this work. The H-H, O-H, and O-O repulsive interaction were not sampled in our training set and were consequently taken from the DFTB+ miomod-hh-0-1 and mio-1-1 parameter sets. All DFTB/ChIMES calculations were subsequently performed with an SP level of theory, using periodic boundary conditions and the same k-point sampling as used for our DFT calculations. The main goal of our DFTB/ChIMES calculations was to determine initial values for diffusion hop activation energies that could then be refined by more computationally expensive DFT calculations. Further details of the ChIMES formalism and our DFTB/ChIMES model are included in the
Supplementary Materials. Refinement of our DFTB/ChIMES model is the subject of future work.
3. Results
We have computed a survey of results from PBE, PBESol, and SCAN by determining optimized ground-state PuO
properties for a variety of spin geometries and with GGA + U (
Table 2). All of our results for all functionals and levels of theory yield lattice constants that are close to the experimental value and lattice volumes that agree within ∼3%. We find that all magnetic states (FM and AF) are essentially isoenergetic for all three functionals, while the diamagnetic state exists over 3 eV per formula unit higher in energy for PBE and PBESol, and over 5.5 eV for SCAN. In addition, the FM (111) state is
eV higher than all AF states for PBE and PBESol, but is
eV lower in energy for SCAN. However, these energy differences are negligible at 298 K (~0.025 eV).
As expected, all three AF calculations with SP yield virtually isoenergetic structures, though there are differences in the lattice constants and in some cases the final lattice volume. In all cases, the AF (1klong) and AF (1ktrans) spin states resulted in a tetragonal lattice distortion of similar magnitude (0.06 Å). In contrast, the AF (3klong) state yields a smaller tetragonal distortion for PBESol (0.02 Å), an orthorhombic unit cell for PBE, and a cubic unit cell for SCAN. In general, PBESol appears to yield the best overall agreement with the experimental lattice constant and volume for the magnetic states studied here, though the differences between PBE and PBESol do not appear to be substantial. We note that can be adjusted in order to yield results close to the experimental lattice volume.
Some additional differences are observed with some of the band gap results. All three AF band gaps from PBESol agree with a value of 2.92 eV. However, both PBE and SCAN yield a slightly higher band gap for AF (3klong) state relative to AF (1klong) and AF (1ktrans), with values of 2.85 vs. 2.80 eV for PBE and 2.75 vs. 2.68 eV for SCAN. PBESol yields a higher gap for the DM state relative to the AF states (3.07 eV), whereas PBE yields a lower gap (2.60 eV), and SCAN yields a metallic system. All three functionals yield a lower gap for FM (111) compared to the AF states, with a value of 2.54 eV from PBESol, 2.60 eV from PBE, and 2.16 eV for SCAN. Overall, the AF results from PBE are generally closest to the experimental value of 2.82 eV, though all three functionals yield reasonable results. Results for AF (1klong) GGA + U calculations from DFTB/ChIMES yielded an accurate lattice constant for a value of eV, albeit for a metallic state. We found that increasing the value of to as high as 12 eV did not create an insulating state. This is possibly due to either the minimal basis set used in DFTB calculations or the omission of greater than two-center terms in the DFTB Hamiltonian. Use of different confining radii and/or basis set functions could improve upon this result, though those effects are beyond our current scope of work.
Finally, in order to test the validity of our SP results, we have performed noncollinear calculations with SOC and the same value of
eV. We find that this set of calculations yields the AF (3klong) SOC state as the energetic minimum, with all three magnetic states yielding a cubic structure and lattice constant very close to that of experiment. However, once again, the AF (3ktrans) SOC and FM (111) SOC states are higher in energy by <0.01 eV, indicating that these spin geometries are essentially isoenergetic. In addition, we compute a FM (111) SOC band gap of 2.50 eV, which is very similar in value to the SP result. However, we compute AF (3klong) SOC and AF (3ktrans) SOC band gaps that are higher than the SP result by 0.17 eV (3.09 vs. 2.92 eV). The SOC relative energy and band gap results reported here differ somewhat from previous PBESol SOC results using the same value for
[
13]. These are possibly due to a combination of the larger k-point mesh used in our study as well as the more stringent convergence criteria used in Ref. [
13].
Based on these results, we choose to proceed with the PBESol functional for the rest of our calculations. PBESol provides accurate equilibrium lattice constants, particularly for the nearly cubic AF (3klong) magnetic state, as well as a band gap only slightly higher than that from experiment. Hence, it appears to provide the best combination of accuracy and computational efficiency, though we note than any of three functionals studied here would appear to provide reasonable results.
Next, we have computed hydrogen point defect formation energies for the int and oct sites using different magnetic states and levels of theory. These were calculated according to the equation:
where the subscript ‘def’ corresponds to either the int or oct site,
is the total int or oct defect system energy,
n is the number of formula units in the PuO
/H system,
is the bulk energy per formula unit, and
is the molecular hydrogen formation energy. The goal here is to assess the effects of these parameters on interactions between the lattice and the hydrogen ion itself. For most spin states, we have fixed our simulation cell at the experimental lattice parameters, which are close the optimized values reported in
Table 2. The only exceptions is the AF (1klong) calculation where the supercell lattice was held fixed at the dimensions reported in the table in order to provide direct comparison to previous H/PuO
collinear SP GGA + U results [
16,
20]. One additional AF (1klong) calculation was performed with the experimental lattice vectors with SP and without GGA + U, which is labeled as AF (1klong-iso).
Our results indicate a relative consistency in the hydrogen interstitial site formation energy,
(
Table 3). For the SP calculations, results without GGA + U yield a fairly narrow spread of
, with a value of 0.63 eV for FM (001) at the low end and 0.71 eV for AF (3klong) at the high end. The effect of using GGA + U appears to be fairly small, with a result of 0.64 eV from AF (3klong) GGA + U. In this case, additional use of SOC has a similar magnitude of effect in the opposite direction, with a result of 0.70 eV for AF (3klong) GGA + U SOC. The tetragonal cell distortion in the AF (1klong) calculation does not appear to have a significant effect on
, though we note that the breaking of the cubic symmetry likely results in the creation of additional interstitial sites and energies. Our DFTB/ChIMES result agrees relatively closely with the spread of results from PBESol with a value of 0.75 eV. Our results at all levels of theory for all magnetic states agree with the previously published SP GGA + U result from PBE of 0.54 eV from Ref. [
16]. The exception presented here is the diamagnetic state, which yields a relatively high
value of 0.97 eV. Our results also show some agreement with the range of PBE results of 0.48-1.19 eV from Ref. [
20], though there is some ambiguity in that publication as some calculations might have been performed on the PuO
unit cell, where system size effects are likely significant.
Given the relative agreement between levels of theory, we choose to continue our study with the more efficient SP calculations, only. Results for the octahedral site formation energy,
, with results from collinear FM (001) and noncollinear AF (3klong) showing agreement within 0.08 eV, with values of 1.25 and 1.33 eV, respectively. Inclusion of GGA + U in the AF (3klong) GGA + U calculation yields a modest increase in
of 0.16 eV to a value of 1.49 eV. We note that our AF (3klong) GGA + U result is smaller than the collinear SP GGA + U results from PBE of 2.24 eV [
16] and 2.12 eV [
20]. These previous PBE results are in relatively good agreement with the DFTB/ChIMES result of 2.45 eV.
Regardless, our results indicate overall that both SOC and GGA + U likely play a relatively small role in hydrogen/PuO lattice interactions. Hence, for the remainder of this work, we choose to perform calculations with SP, only and omit these two levels of theory for the sake of computational efficiency. This allows us to more efficiently explore the kinetics of hydrogen diffusion pathways as a function of magnetic state.
We now investigate the kinetics of three fundamental pathways in the H/PuO
system (
Figure 2), namely:
Diffusion from one int site to its nearest neighbor int site (where an O–H bond is broken), which we label int-int.
Diffusion between int sites containing an O–H bond rotation (where the bond remains intact), which we label int-int2.
Diffusion hops from the oct site to the closest int site, which we label oct-int.
We have not computed int-oct diffusion (from an interstitial to an octahedral site) due to the large increase in energy for octahedral absorption, indicating that this pathway is likely only at extreme temperatures. In addition, oct-oct diffusion is highly unlikely to occur, given the closer proximity of low energy int sites.
We now approximate the activation energies,
, for these pathways using a two-step approach. First, we have computed minimum energy pathways (MEP) from our DFTB/ChIMES model using the climbing image nudged elastic band (NEB) method [
43]. All calculations were run within the Atomic Simulation Environment python interface [
44,
45], with the image dependent pair potential method [
46] used to generate an initial guess for each pathway. Here, our SCC convergence was held to
eV (
au) and the NEB forces were converged to
eV/Å. In addition, the simulation cell was held fixed at the experimental PuO
density with cubic symmetry. Decreasing the SCC convergence to
eV (
au) and the NEB force convergence to
eV/Å resulted in an increase in the activation energy of 0.01 eV for the int-int2 rotation. Next, we performed constrained geometry DFT optimizations starting from the DFTB/ChIMES computed transition state geometries using PBESol with different spin states, where the hydrogen ion was held fixed and the surrounding PuO
ionic positions were allowed to relax. In doing so, we were able to approximate
for several different spin states relatively efficiently at a DFT level. All PBESol optimizations were kept at the experimental PuO
cubic lattice, excluding AF (1klong) where the ionic coordinates and lattice vectors were scaled to the optimized tetrahedral distorted cell in order to facilitate comparison with previous DFT NEB results [
20]. More exact calculations of MEPs in H/PuO
systems is the subject of future work.
Our results indicate a high degree of consistency between results for different spin states (
Table 4). All magnetic states from PBESol conducted on a cubic simulation cell yield an
of ~1.1 eV for the int-int transition, while the diamagnetic calculation yields a value of only 0.36 eV. DFTB/ChIMES yields a value in agreement with the PBESol magnetic results, with a value of 1.03 eV. Additional constrained optimization calculations were performed with PBE in order to examine the effect of choice of functional on our results. We find that results from both the diamagnetic and AF (3klong) states yield close agreement with results from PBESol, with values of 0.32 eV and 1.10 eV, respectively. Our result for the tetragonally deformed AF (1klong) geometry computed with PBE yields a value of 0.15 eV along the short lattice length, very close to the previously reported PBE activation energy of 0.13 eV for this transition [
20].
Given the relative agreement between results from different magnetic states with a cubic simulation cell, we continue with the AF (3klong) spin geometry for PBESol, only. We find that for the int-int2 transition DFTB/ChIMES determines a relatively large relative to PBESol, with values of 0.36 eV vs. 0.06 eV. In contrast, we observe improved agreement for oct-int, where DFTB/ChIMES yields an activation energy of 0.27 eV compared to the value of 0.36 eV from PBESol. It is possible that the agreement of DFTB/ChIMES with PBESol could be improved with additional and more accurately computed training data, which is the subject of future work. We choose to use the PBESol AF (3klong) estimates for for the remainder of this work due to their likely higher accuracy than those from DFTB/ChIMES.
We estimate quantum effects on hydrogen diffusion by treating the diffusion pathways as a one-dimensional quantum double-well with a parabolic barrier, using an approach adapted from Ref. [
47]. We start with the following piece-wise quadratic function:
Here, k is the spring constant for each well, is the spring constant for the barrier, is the position of the reactant minimum, is the position of the product minimum, is the parabolic barrier height, is the offset of the product well, and and define the reactant (left) and product (right) regions. For our work, we have chosen in order to ensure a smooth potential throughout, where . The parameter is set to the approximate value from the AF (3klong) calculations. The potential is centered at zero for all three transitions studied here. For the symmetric int-int and int-int2 transitions, , and , where is the total reaction path distance as determined from our DFTB/ChIMES NEB calculations. Furthermore, we compute the reaction/product positions as and . For the asymmetric oct-int transition, the same equations apply with the exception of and .
We can estimate the quantum tunneling rates by first assuming a quantum harmonic oscillator (Hermite polynomial) basis set and solving for the energy eigenstates. Care was taken to ensure convergence for each transition in terms of basis set size as well as the numerical limits of integration. Results for the int-int transition indicate small splittings between energy levels (i.e., <1 meV) (
Figure 3), typical for a symmetric double well with a relatively large barrier where the particle tends to remain in a single well and tunneling transmission coefficients are small. The probability of hydrogen being found between wells becomes more pronounced as the eigenstates approach the top of the potential barrier, with the state just below the barrier yielding a splitting of slightly less than 0.01 eV. In addition, we compute a zero-point energy correction in each potential well of ~0.06 eV, and that the energy levels below the parabolic barrier tend to be spaced by ~0.13 eV.
The tunneling transmission coefficient for each eigenstate can be computed using the analytical expression for a parabolic barrier [
48]:
The thermal reaction rate for hydrogen diffusion for a given pathway is then computed from the following expression:
The summation is over the
n computed eigenvalues
and eigenstates
of the system, and
is the probability of finding the particle on the reactant side, namely:
Similar to Ref [
47], we choose to define the reactant region limit for a given eigenstate,
, as its classical turning point. We note that the computed thermal reaction rate will have some dependence on definition of
, though the classical turning point is physically intuitive and consistent with well-known semi-classical approaches such as the Wentzel-Kramers-Brillouin (WKB) approximation.
Additional accuracy for our model double-well potentials could involve determination of flux-flux autocorrelation functions [
49] and/or time-dependent wave-packet propagation to determine the scattering matrix [
47,
50]. Full determination of quantum nuclear vibrational effects in the condensed phase would require simulations using path integrals (e.g., Ref. [
51]), which would necessitate the development of an accurate interatomic potential due to the extreme computational expense of these simulations [
52]. Regardless, results from Ref. [
47] indicate that the eigenstate transmission approach we have adopted here likely yields a high degree of accuracy for proton double-well transmission barriers as low as ~0.3 eV with respect to scattering matrix calculations.
We now make comparison to classical rate constants estimated from our approximate
values from the PBESol AF (3klong) calculations. Classical diffusion rates were estimated from a modified version of the Eyring equation:
Here, we have assume that the free energy of activation is equal to our estimate for
and that the Eyring transmission coefficient is equal to one. This approach yields a kinetic prefactor of
at 300 K, which is of the same magnitude as DFT calculations for hydrogen ion diffusion on transition metal oxide surfaces (e.g., Ref. [
53]).
Our results for the diffusion rates (
Figure 4) indicate that int-int diffusion is likely low even at 300 K, where we compute a
value of
and a
value of
. We predict that the rate increases by nine orders of magnitude at 600 K, with values of
and
. These increase by an additional four orders of magnitude at 1000 K, with a values of
and
.
In contrast, the oct-int transition is relatively frequent at 300 K due to its relatively low barrier, with values of and . These values increase more modestly with increasing temperature, with values of and at 600 K and at and 1000 K. Nonetheless, hydrogen absorption into an octrahedral defect site is unlikely at the conditions studied here, given its high formation energy of eV.
Finally, we compute relatively rapid int-int2 O–H bond rotations for all temperatures studied here, with values of and at 300 K. The temperature effect on the rate is expected to be small, given the relatively low barrier of 0.06 eV (~700 K). In this case, we compute rates of and at 600 K, and and at 1000 K. This transition is likely to only play a secondary role in hydrogen mass transport through the PuO lattice, given that the O–H bond stays intact throughout this pathway.
The ratio of the quantum to classical rate indicate significant quantum effects for the int-int and oct-int transitions (
Figure 5), both of which involve the breaking and/or forming of a covalent O–H bond. The oct-int diffusion pathways appears to show the strongest enhancement at low temperature due to quantum tunneling with a
ratio of ~180, compared to the int-int ratio of ~19. However, our approach uses an identical prefactor for all pathways at all temperatures, whereas harmonic transition state theory could indicate that the prefactor for the int-int path is considerably higher than that of oct-int due to the presence of a high frequency O–H bond in the int absorption site. Regardless, our results provide reasonable estimates for the persistence of quantum nuclear vibrational effects even at high temperature. Both int-int and oct-int pathways show a quantum enhancement of approximately a factor of three at 1000 K. The relatively low
of the int-int2 transition yields more modest quantum effects at 300 K of slightly less than
enhancement, which approaches a value close to one at 1000 K. Overall, our results indicate that quantum nuclear vibrational tunneling effects could enhance hydrogen diffusion in plutonium oxide materials by more than one order of magnitude at ambient conditions. Hence, it is likely that these effects need to be included in order to make accurate predictions for experiments. More accurate quantification of these effects can accessed through the development of reactive molecular dynamics models, which is the subject of future effort.