# Estimates of Quantum Tunneling Effects for Hydrogen Diffusion in PuO\({}_{2}\)
## 2 Computational Details  
PuO\({}_{2}\) has been determined experimentally to have a cubic lattice with \(Fm\bar{3}m\) 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\({}_{2}\) 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 \(6\times 6\times 6\), energy convergence criteria of \(10^{-6}\) eV and force convergence criteria of \(10^{-2}\) eV/A, which yielded system pressures far below 1 GPa. Tests with PBE using a planewave cutoff of 600 eV, k-point mesh of \(10\times 10\times 10\), energy criteria of \(10^{-8}\) eV, and force criteria of \(10^{-3}\) eV/A indicated convergence of energies within \(\sim 3\times 10^{-3}\) eV per formula unit and lattice constants within 0.01 A. 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^{-6}\)A. 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 \(2\times 2\times 2\), energy criteria of \(10^{-5}-10^{-6}\) eV, and force criteria of \(5\times 10^{-2}\) eV/A. Use of a linear response approach [35] to compute \(U_{\text{eff}}=U-J\) 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 \(U_{\text{eff}}=6.00\) eV for our calculations with \(U=6.00\) eV and \(J=0.00\) 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].