A first-principles study of hydrogen surface coverage on δ-Pu (100), (111), and (110) surfaces

Hydriding corrosion of plutonium leads to surface cracking, pitting, and ultimately structural failure. Laboratory experiments demonstrate that hydriding begins on the surface or near subsurface of plutonium. However there has not yet been a systematic evaluation of hydrogen surface coverage on plutonium. In this work, we compute the surface energies of the low facet surfaces of face-centered cubic δ-Pu. The adsorption free energies of expected hydrogen structures at low and high coverage are presented, along with the likely progression for filling sites as the H2 partial pressure increases. Implications for near-equilibrium pressure hydride nucleation and non-equilibrium millibar pressure hydriding are discussed.


I. INTRODUCTION
Plutonium's unique material properties make it useful as a heat source in power plants and in radioisotope thermoelectric generators. 1 However, plutonium is highly reactive and will corrode if exposed to atmospheric gases, even in trace amounts. Hydriding in face-centered cubic (fcc) δ-Pu occurs if the hydrogen concentration exceeds a temperature-dependent solubility limit. As it incorporates hydrogen to form plutonium dihydride (PuH 2 ), the plutonium lattice expands by approximately 60 volume %, causing cracks and fissures in the plutonium surface. Moreover, PuH 2 is pyrophoric, catalyzing rapid, violent exothermic oxidation. As a result, uncovering the mechanisms by which plutonium corrodes can improve safety regulations and eliminate toxic waste.
Conducting experiments with plutonium is difficult due to its toxicity and radioactivity.
In addition, the complex phase diagram of plutonium makes it difficult to grow the single crystal materials necessary to characterize surface properties such as the surface energy. In particular, hydriding experiments are difficult since a PuO 2 layer of nonuniform thickness is present on plutonium surfaces that have been exposed to air, providing variable protection to the underlying plutonium and introducing an element of randomness into the observed hydriding induction times. 2 Lastly, experiments alone cannot provide a detailed atomicresolution hydriding mechanism nor identify ultimately how it can be controlled.
Quantum simulation methods, such as density functional theory (DFT), are necessary to model the bond breaking inherent in hydrogen dissociation and adsorption. [3][4][5][6] DFT remains one of the most popular modeling methods in condensed matter physics, computational chemistry, and materials science, and has been widely used in studies of the physical and structural properties of δ-Pu (e.g., Refs. 7,8). While a comprehensive explanation of plutonium's magnetism has yet to be offered, plutonium does not exhibit a bulk magnetic moment 9 and anti-ferromagnetic (AFM) models yield structural properties comparable to experiments. 10 Surface hydriding studies are generally framed in terms of low (single molecule per surface unit) vs. high (nearing full monolayer) coverage. Low coverage conditions will be expected at hydrogen partial pressures P H 2 well below the equilibrium pressure P eq H 2 , P H 2 P eq H 2 . Previous DFT studies of adsorption of a single H atom or H 2 molecule on clean (100) and (111) surfaces of δ-Pu have all been conducted at low coverage conditions. [3][4][5][6][11][12][13] High surface coverage will be expected when P H 2 P eq H 2 . The equilibrium pressure P eq H 2 at room temperature, has not been experimentally verified but was recently predicted from simulations to be O(10 −19 ) bar. 14 Even very small pressures (e.g. 1 mbar) will exceed P eq H 2 by many orders of magnitude resulting in rapid adsorption and, unless the surface has a protective layer of PuO 2 , barrierless hydriding. Immediate hydriding at millibar pressures has been observed by experiments. 2,[15][16][17][18] We should therefore anticipate that under real world conditions, adsorption onto the surface from the gas phase may occur more quickly than surface adatoms can burrow from the surface into the bulk. However, to date there have been few systematic studies of Pu surface hydriding including determination of initial mechanisms as well as studies of middle-ranged surface coverage and how hydrogen can accumulate on a given facet.
In the present work, we use DFT calculations to determine hydrogen adsorption energies ranging from low to high surface coverage on low-index facets on δ-Pu. We first use a simple two-state model to discuss relevant hydrogen surfaces coverages in terms of adsorption energetics when P H 2 ≈ P eq H 2 . This simple model helps place bounds on the numbers and types of adsorptions and their resulting energies that are likely relevant to initiating surface hydriding on different crystalline facets. We quantify bulk absorption of hydrogen in δ-Pu and PuH 2 at varying hydrogen concentrations, including the formation of point defects in PuH 2 .
We then use DFT calculations to evaluate the surface energies of the (100), (111), and (110) surfaces to establish the relative abundance of each surface in polycrystalline materials. We explore three AFM geometries for each crystalline facet by alternating magnetic spins along (100), (111), and (110) planes. These results are also discussed in terms of the thermodynamic stability and morphologies of different nanocrystals. 19 For the lowest energy AFM configuration of each surface, we next compute the adsorption energies of anticipated hydrogen structures with increasing surface coverage and compare to bulk absorption energies.
Our efforts help determine an understanding of how increasing coverage affects subsequent hydrogen adsorption on the various surfaces of δ-Pu and will provide insight into possible hydriding mechanisms observed in the bulk.

II. COMPUTATIONAL DETAILS
Within DFT, f -electron correlations in bulk plutonium and plutonium hydride can be accounted for using either an on-site Coulomb repulsion (e.g. GGA+U) or orbital polarization, while the magnetic state can be computed through spin polarization (SP) or spin-orbit coupling (SOC) calculations. [20][21][22][23][24][25][26][27][28][29][30][31] While spin-orbit coupling (SOC) is essential to compute material properties of PuO 2 , 25,32,33 it is less critical in simulating H-Pu interactions in either the surface or bulk. 25,29 SP computed energies of low-index plutonium surfaces differ from SOC computed surface energies by only 7%. 8 In addition, adsorption energies of isolated hydrogen atoms on the (100) and (111) surfaces of plutonium differ by only −0.03 eV on average when computed with SP vs. SOC. 6 As a result, we choose to employ the less computationally intensive GGA+U and SP methods in this work.
All DFT calculations discussed here were performed with the VASP code 34-36 using projector-augmented wave function (PAW) pseudopotentials 37,38 and the Perdew, Burke, and Ernzerhof exchange correlation functional (PBE). 39 The cutoff energy for the planewave expansion was set to 500 eV and the wave function was converged to within 10 −5 eV.
Partial occupancies of the electronic states were set with a fourth-order Methfessel-Paxton smearing 40 with a width of 0.3 eV.
Bulk δ-Pu and PuH 2 were modeled using 32 plutonium atoms and periodic boundary conditions. The magnetic moments of Pu atoms were treated using collinear SP and were alternated along (100) planes to give antiferromagnetic (AFM) configurations. A 5 × 5 × 5 Monkhorst-Pack mesh 41 was used for all bulk calculations, whereas we used a mesh of 5×5×1 Slabs of δ-Pu were created with exposed (100), (111), and (110) surfaces. Each slab consists of 4 layers of plutonium atoms with 30 Å of vacuum. Atoms in the bottom two layers were held fixed in place, while atoms in the top two layers were allowed to relax.

III. TWO-STATE MODEL
We consider surface coverge under equilibrium conditions using a simplified two-state model. For a system of N H hydrogen atoms partitioned between N bulk bulk sites and N surf surface sites, the surface coverage is where N surf H is the number of hydrogen atoms adsorbed to surface sites. The absorption energy at each bulk (surface) site is E bulk (E surf ), and the difference between these energies ignoring constant terms.
At constant temperature T , the average number of atoms adsorbed on surface sites is where Ω is the microcanonical density of states, Q is the canonical partition function, and the dependence of Ω and Q on N H , N bulk , and N surf have been omitted for clarity. The Here, the second probability related to N bulk H was computed via Sterling's approximation for the sake of simplicity.
We consider slabs with N surf = 100 sites and N bulk ranging from 10 4 − 10 7 . The number of hydrogen atoms is set to N H = 10 −3 N bulk to match the upper bound for H absorption in δ-Pu at room temperature. 14 At each system size, we then vary the value of ∆E from highly positive values (corresponding to strongly favorable bulk absorption) to highly negative (corresponding to strongly favored surface adsorption) in order to determine the range of ∆E that will permit high surface coverage.  This simplified two-state model of adsorption illustrates convenient guidelines to evaluate surface coverage at equilibrium. Sites with adsorption energies within ±2 k B T of the bulk absorption energy will be populated at approximately the bulk ratio, N H /N P u . Surface sites that are 10 k B T higher in energy than bulk adsorption will have negligible population. Sites that are 10 k B T lower in energy will be entirely filled. We can thus use ±10 k B T as an approximate bound for gauging how absorbed hydrogen atoms partition between surface and bulk sites. The first step in applying this framework to H surface coverage on δ-Pu is to compute bulk absorption energies, which is where we now turn.

IV. RESULTS
A. Hydrogen absorption in bulk δ-Pu and PuH 2 The energy released as hydrogen atoms are absorbed into bulk δ-Pu from the gas phase where e bulk P u is the energy of a plutonium atom in a defect-free fcc lattice, e H is the energy of a hydrogen atom in a gas-phase hydrogen molecule, and N P u and N H are the number of plutonium and hydrogen atoms in configuration i, respectively. The absorption energy reveals that it is energetically favorable for an isolated hydrogen atom in bulk δ-Pu to absorb at an O site (∆E bulk O1 = −2.92 eV) rather than at a T site (∆E bulk T 1 = −2.71 eV). The difference ∆E bulk T 1 − ∆E bulk O1 = 0.21 eV (8 k B T at room temperature) and the corresponding Boltzmann factor, exp (-8), shows that the fraction of absorbed hydrogen occupying T sites is O(10 −4 ).
There are three possible nearest neighbor N H = 2 clusters in bulk δ-Pu: neighboring occupied O sites (O2), neighboring occupied T sites (T2), and an occupied O site next to an occupied T site (T1O1). The O2 absorption energy ∆E bulk O2 is 0.25 eV (10 k B T ) greater than the energy of two independent O1 absorbates, 2∆E bulk O1 , so the equilibrium population of O2 configurations will be O(10 −5 ) of the absorbed hydrogen. The T1O1 and T2 absorption energies are even higher, 21 k B T and 27 k B T above the O1 absorption energy, respectively, indicating that the equilibrium population of either is negligible.
The preceding analysis demonstrates that ∆E bulk O1 is a convenient reference for evaluating absorption energies. Normalizing ∆E bulk i by N H and subtracting the reference absorption ∆∆e bulk i will be negative for configurations that are more stable than a set of N H atoms distributed across isolated O sites, and positive for configurations that are less stable.
is 3.24 eV. Similarly, the energy released by an O site interstitial where E s P u is the DFT energy of a plutonium slab with two exposed s-oriented surfaces and A s is the exposed surface area. Surface energies are reported in Table I  We find that in general the AFM configurations have lower surface energies than the corresponding FM surfaces and that the surface energy is minimized by the AFM configuration with planes of magnetization that are parallel to the surface. As a result, we would

C. (100) Surface monolayer formation
On a surface of δ-Pu, the adsorption energy is where E s i is the DFT energy of a plutonium slab with s-oriented surfaces and with configuration i of surface hydrogens. ∆∆e s i is defined similarly to Eq. 5 The (100)  We now sequentially fill the H and B sites on the (100) surface and compute ∆∆e

D. (111) Surface monolayer formation
On the (111) surface of δ-Pu, hydrogen atoms can adsorb at fcc or hexagonal close-packed (hcp) sites, which differ only in the positioning of plutonium atoms in the first subsurface layer. If a hydrogen atom moves from a surface hcp site into the bulk, it will occupy a T site. If a hydrogen atom in a surface fcc site moves into the bulk, it will occupy an O site.
The (111) slab used here has 8 surface Pu atoms, and therefore has 8 fcc and 8 hcp sites.
We find that the adsorption energy of a single H atom at either site are indistinguishable,  lowest adsorbtion occurs at fcc sites (see Fig. 8N) and is only a few k B T more favorable than bulk adsorption, ∆∆e As surface coverage on (110) increases, we consider surfaces that have all or all but one of a single site type filled. In this limit, LB sites are the only ones that still have more favorable adsorption than the bulk. If all but one LB site are filled (see Fig. 8O), ∆∆e  After the LB sites are filled, hydrogens occupy fcc sites and the adsorption energy increases more shallowly, but still approximately linearly.

V. DISCUSSION
A picture of near-equilibrium hydriding begins to emerge as we examine adsorption energetics on all three surfaces. On the (100) surface, the lowest energy sites are H sites, which are 1-2 k B T /H higher in energy than bulk O sites. Accordingly, the two-state model results suggest that the fraction of populated H sites will be comparable to the bulk absorption, O(10 −3 ). The (111) surface will have no substantial surface coverage since the lowest energy structures are all > 7 k B T /H higher in energy. On the (110) surface, adsorption in isolated LB sites is −15 k B T /H relative to the bulk, suggesting that some LB sites will always be filled. Even with the majority of LB sites filled the adsorption energy increases to only −3 k B T /H. The two-state model, in turn, shows that surface sites with −3 k B T adsorption energies will be only be 10% filled, though this picture is somewhat oversimplified. Regardless, LB sites, though few in number due to the high surface energy of the (110) surface, could serve as nucleation sites for equilibrium hydriding due to their high stability once occupied.
We can also develop a different picture for non-equilibirum hydriding at millibar H 2 pressures by comparing our adsorption energy studies to previously published values of H 2 molecular adsorption and dissociation and H ion subsurface penetration on the (100) and (111) surfaces. Published data of this type is not available for the (110) surface and hence we limit our discussion to these two surfaces. On (100), an adsorbed H 2 molecule is held to the surface only weakly, 0.15 eV (6 k B T ). 3 Computed dissociation energy barriers vary, 0.22 eV (9 k B T ) 6 to 0.49 eV (19 k B T ) 13 to 0.78 eV (30 k B T ), 3 but are all higher than the desorption barrier, suggesting that hydrogen molecules desorb and re-adsorb many times before dissociating. Once dissociated, the barrier to recombine is 1.44 eV (56 k B T ), while the energy barrier for a surface hydrogen atom to burrow into a subsurface T site is 0.93 eV (36 k B T ). 13 Comparing the barriers for burrowing and recombination suggests, in turn, that surface hydrogen atoms are more likely to move into the bulk than recombine into an H 2 molecule.
A similar but distinct picture emerges for the (111) surface. For an adsorbed H 2 molecule to detach from the surface requires 0.13 eV (5 k B T ). 4 Dissociation may require less energy, 0.06 eV (2 k B T ) 13 , or comparable energy, to 0.12 eV (5 k B T ) 6 to 0.31 eV (12 k B T ). 4 Either way, both desorption and dissociation can be considered barrierless. Once dissociated, the recombination barrier is 1.36 eV (53 k B T ). 13 For comparison, the energy barrier for a surface hydrogen atom to burrow into the subsurface layer is half, 0.66 eV (26 k B T ). 13 Moreover, the high adsorption energies computed in this work provide a high driving force for diffusion into the bulk, which is possible from both fcc and hcp sites. Overall, the (111) is a likely candidate for non-equilibrium hydriding since dissociation is barrierless and there are many pathways to subsurface diffusion.

VI. CONCLUSION
We have evaluated the surface energies of the low index facets of δ-Pu. We find that the lowest energy surfaces have anti-ferromagnetic geometries with magnetic spins alternating along planes parallel to the surface. Of these, the (111) surface has the lowest surface energy and so would be the dominant surface of a macroscopic single crystal material.
Polycrystalline materials will exhibit multiple surfaces; we predict the (111) will be the most abundant and the (110) the least abundant.
We have also evaluated hydrogen adsorption energies ranging from low to high surface coverage on each surface. We find that the (111) surface is the least attractive to adsorbed hydrogen, while the (110) surface is the most attractive. This suggests two mechanism for hydriding: equilibrium hydriding nucleated by the most attractive LB sites and non-equilibrium hydriding facilitated by barrierless dissociation and high driving forces for subsurface diffusion on the (111) surface.
Future model development of plutonium surfaces needs to account for restructuring of surface atoms and the presence of surface defects present in real plutonium. Additionally, subsurface absorption energies as a function of depth below the surface, particularly the number of layers required before absorption energies are equal to bulk values, will further our understanding of the role that surface and subsurface effects play in near-equilibrium hydriding. Overall, our results yields a relatively simple description of the onset of hydriding in δ-Pu that can be used to provide constraints for macro-scale models of material aging and degradation. Our efforts can be extended to include interactions with plutonium oxide surfaces, which could yield advances in the way Pu-containing devices are designed in future applications.

SUPPLEMENTARY MATERIAL
See supplementary material for the following information: