Uranium (U) is an actinide element exhibiting delocalized f-electrons that exists in three solid phases: α (face-centered orthorhombic), β (tetragonal with a 30 atom unit cell) and γ (body-centered cubic) (Yoo et al., 1998). Below 935 K, U exists in the α phase (Söderlind, 1998). Although uranium is alloyed with Mo or Zr in order to stabilize the preferred body-centered cubic phase at lower temperatures, most notably for nuclear fuel applications, the α phase persists either in decomposed phase regions (as in the case of U-Mo fuel) or in the low-temperature periphery of the fuel slug (as in the case of U-Zr type fuel). Because of its unique crystal structure, α-U exhibits anisotropic material properties, for example, thermal expansion (Lloyd and Barrett, 1966) and elastic constants (Fisher and McSkimin, 1958). One relevant behavior is the anisotropic irradiation growth of α-U (Leggett et al., 1963), which can lead to the tearing along grain boundaries and the formation of crystallographically aligned pores. These porosities are believed to form due to the interaction of lattice defects and stresses induced by irradiation and a temperature gradient and result in significant volume changes of α-U (Paine and Kittel, 1958; Leggett et al., 1963). A full understanding of the relationship between point defects, lattice strains, and the generation of pores, has not been established either experimentally or via computational methods. However, a number of computational studies have been performed into fundamental properties of uranium (Jokisaari, 2020).

Density functional theory (DFT) is an essential part of computational materials science, addressing a variety of problems in materials design and processing on a fundamental level. Several examinations of U via DFT have been performed on its orthorhombic and body-centered cubic (bcc) structures. (Taylor 2008) used a projector augmented wave (PAW) pseudo-potential to calculate the lattice constants of α-U and γ-U along with the bulk modulus of both phases and the surface energy of a single surface in α-U. (Xiang et al., 2008) also utilized a PAW pseudo-potential to perform an analysis of lattice constants in the α and γ phases, as well as an analysis of defects in γ-U. Beeler, et al. calculated the lattice constants and elastic constants of α, β, and γ-U (Beeler et al., 2013), in addition to the point defect properties in both α and γ-U (Beeler et al., 2010). (Huang and Wirth, 2011; Huang and Wirth ,2012) calculated intrinsic and extrinsic point defect formation energies and migration barriers in α-U. These DFT investigations showed a general, if imperfect, agreement with the experimental bulk moduli, lattice and elastic constants, and vacancy formation energy (Barrett et al., 1963; Matter et al., 1980; Yoo et al., 1998).

Key areas lacking from the DFT-based literature on α-U include an analysis of temperature-dependent bulk properties, including thermal expansion and heat capacity, and temperature-dependent defect properties. Extrapolation of 0 K properties to compare to elevated temperature experiments can be ill-advized, particularly for complex crystal structures or for point defects. However, ab initio molecular dynamics (AIMD) allows for quantum mechanical-based calculations to be performed at non-zero temperatures. AIMD has been utilized to study a variety of systems including liquid phase diffusion in Al-Si (Manga and Poirier, 2018), the adsorption energy of Fe on TiN surfaces (Wang et al., 2010), NaCl dissolution in water (Timko et al., 2010) and finite temperature phonon dispersion curves in bcc Zr and bcc (Li Hellman et al., 2011). Hood et al. (2008) have previously utilized AIMD to study the equation of state of U and the variation in density of states for the liquid phase of U at two unique temperatures. Hood, et al. utilized a unique pseudo-potential that was presented in that same manuscript (Hood et al., 2008). Söderlind et al. (2012) utilized self-consistent ab initio lattice dynamics (SCAILD) to study the high-temperature stabilization of the γ-U phase by calculating phonon modes at 1,100 K. There have been several investigations of the free energy via the temperature-dependent effective potential (TDEP) technique (Hellman et al. 2013; Bouchet and Bottin, 2017; Kruglov et al., 2019; Castellano et al., 2020; Ladygin et al., 2020), and one recent AIMD investigation into the defect formation energies in γ-U at high temperatures (Beeler et al., 2020). Despite numerous investigations into metallic uranium, AIMD has not been deployed to study temperature-dependent properties of α-U, such as thermal expansion or defect formation energies.

In this work, AIMD simulations are performed to calculate the equilibrium volume of α-U as a function of temperature. Utilizing these equilibrium volumes, the three lattice constants of α-U are determined, yielding anisotropic thermal expansion calculations. Additionally, point defect formation energies for interstitials and vacancies are calculated as a function of temperature. The lattice strain from the introduction of these point defects is calculated, providing insight into the fundamental stages of irradiation growth. This is the first AIMD investigation of point defects in α-U.

Systems are investigated using the Vienna ab initio Simulation Package (VASP) (Kresse and Hafner, 1993; Kresse and Hafner, 1994; (Kresse and Fürthmuller, 1996a; Kresse and Fürthmuller, 1996b). The projector augmented wave (PAW) method (Blöchl, 1994; Kresse and Joubert, 1999) is utilized within the density functional theory (Hohenberg and Kohn, 1964; Kohn and Sham 1965) framework. Calculations are performed using the Perdew-Burke-Ernzerhof (PBE) (Perdew et al., 1996; Perdew et al., 1997) generalized gradient approximation (GGA) density functional implementation for the description of the exchange-correlation. A uranium PAW pseudo-potential with the 6s²6p⁶5f³6d¹7s² valence electronic configuration and a core represented by [Xe, 5d, 4f] is utilized. Methfessel and Paxton’s smearing method (Methfessel and Paxton 1989) of the first order is used with a width of 0.1 eV to determine the partial occupancies for each wavefunction. A (Monkhorst and Pack, 1976) 1 × 1 × 1 k-point mesh was utilized for Brillouin zone sampling. Uranium is assumed to be non-magnetic, in accordance with both experiments (Lide, 2000) and previous simulations (Beeler et al., 2013), and as such calculations are non-spin-polarized. The Hubbard U term is not applied to the f electrons, in accordance with previous computational studies that utilized DFT on α-U (Taylor, 2008; Huang and Wirth, 2011; Huang and Wirth, 2012; Beeler et al., 2013). The precision is set to normal and the energy cutoff is increased to 300 eV. The electronic self-consistent loop exit criterion is set to 10⁻⁴ and the precision for projectors in real space is increased to −10⁻⁴. A 180 atom supercell (5 × 3 × 3 unit cells) is utilized for all simulations. Such a supercell is sufficiently large to preclude defect-defect interactions (Huang and Wirth, 2011). Dynamics were carried out in the NpT ensemble with Parinello-Rahman dynamics and a Langevin thermostat to control temperature. The temperature range investigated is from 100 to 800 K, in increments of 100 K, to span nearly the entire range of stability of the α phase of U. The Langevin friction coefficient for both the atomic and lattice degrees-of-freedom was set to 5 ps⁻¹, with a fictitious mass for the lattice degrees of freedom of 500 amu. The timestep is set to 2.0 fs, and the initial structures are equilibrated at the target temperature for 8 ps. The trajectories of the total supercell energy, pressure, and individual lattice vectors are analyzed as a function of timestep to ensure that an equilibrated state has been reached and each of these quantities is oscillating around a given value, with the pressure oscillating around zero. Subsequently, the simulations are further equilibrated for 2000 timesteps (4 ps) for the determination of equilibrium properties. In order to obtain average properties over the ab initio molecular dynamics simulation, the energies and lattice constants for the final 1500 timesteps (3 ps) are extracted and averaged. To ensure statistical significance of the results, ten unique simulations are performed for each simulation setup. An example trajectory of the individual supercell lattice constants as a function of time at 300 K is shown in Figure 1. The standard deviation in the a, b, and c equilibrium lattice constants from this single simulation are 0.01, 0.02, and 0.01 Å, respectively. It should be emphasized that this was a single simulation, and ten such simulations were performed to obtain the averaged equilibrium structures and energies. Error bars are excluded from figures in the results that pertain to lattice constants, due to the small deviations observed in the dataset.

The thermal expansion is calculated both for individual lattice constants, and as a volume averaged value, via Eq. 1: α = L − L 0 L 0 (1)

where α is the thermal expansion, and L is the length of the individual lattice constant or averaged lattice constant (cube root of the volume) at a given temperature, and L₀ is the reference length at a given temperature. The formation energy of point defects is calculated via Eq. 2: E f = E * − n ± 1 n × E 0 (2)

where E^* is the energy of a system with a defect (with n ± atoms), n is the number of atoms in the defect-free system and E⁰ is the energy of a defect-free system with n atoms. Eq. 2 utilizes (n + 1) for interstitials and (n−1) for vacancies. The molar heat capacity is determined from Eq. 3: C P = ( δ H δ T ) P (3)

where H is the enthalpy (potential plus kinetic energy), T is the temperature, and it is emphasized that this is determined for a constant pressure, which in these simulations is zero pressure. The total heat capacity is a sum of the lattice (Eq. 3) and the electronic components of the heat capacity. The electronic heat capacity coefficient ( γ e , where the electronic heat capacity is γ e × T) has been experimentally reported to be 10.12 mJ/mol-K² (Marchidan and Ciopec, 1976; Schachinger and Lamprecht, 1989), and this value is utilized here to determine the total heat capacity as a function of temperature.

This work indicates several lattice and volume effects that are strongly dependent upon the temperature of the system. Excepting the heat capacity, all of these calculated properties seem to indicate two regions of behavior, one below 400 K and one above 400 K. Also, it should be re-emphasized that the heat capacity utilized an experimental data point for the determination of the electronic heat capacity. Thus, from the purely computationally obtained data, there are differences in lattice constant, defect-induced lattice strain, and defect formation energy, in these two temperature regimes. To briefly summarize, thermal expansion/contraction in the a/b direction occurs rapidly from 100 up to 400 K, after which a more linear and gradual change in the lattice constant takes place. Defect formation energies decrease in magnitude up to 400 K, above which they tend to slightly increase. Defect lattice strains for interstitials decrease/increase in the a/b direction up to 400 K, above which they maintain essentially the same magnitude. A similar trend in defect lattice strains for vacancies is present, but the magnitude is significantly less and potentially swamped by thermal fluctuations, and thus the statistical certainty of the finding is reduced.

In an attempt to identify underlying changes in the atomic and electronic structure as a function of temperature, two additional analyses were performed. The first was a coordination analysis, in which the radial distribution function (rdf) was constructed for a snapshot of the equilibrated α-U structure at 100, 400, and 800 K. If a modification in the bond lengths is occurring that can yield differences in structure, such a modification would present itself as peak shifts in the rdf. This coordination analysis is presented in Figure 8, but yielded no discernible trends, and the general bond length and peak broadening confirm what is observed in Figures 2, 3 outlining thermal expansion. Specifically, there is general peak broadening with temperature, and the second nearest neighbor peak (in the b direction) slightly shifts toward a smaller distance. The second analysis involved the electronic density of states (DOS). If the character of valence electrons participating in the bonding is changing as a function of temperature, the electronic DOS would identify such changes, and the electronic DOS is known to slightly change as a function of the phase of uranium Beeler et al. (2013). The electronic DOS was calculated, and only minimal variation in the electronic structure was observed as a function of temperature. There was no indication that two regions of electronic structure behavior exist for the α phase. Due to the lack of identifiable structure variation as a function of temperature, likely due to statistical noise in the sampling. These results are not displayed here, as the DOS can be obtained from the literature Hood et al. (2008), Beeler et al. (2013). The underlying cause of the apparent two-region behavior is unclear and warrants further computational investigation, in addition to experimental studies on high-purity uranium to either confirm or refute the findings in this work. It should be emphasized that these results are statistically significant, but are still susceptible to inherent errors in the modeling methodology, such as the inability for DFT to predict the lattice and elastic constants of α-U with sufficient accuracy at 0 K.

This work has identified that individual point defects can produce somewhat significant lattice strain, but it is unknown if such point defects could be responsible for experimentally observed macroscopic shape change. Comparison of the calculated direction of the defect-induced strains to the experimental literature on irradiation growth provides contradictory information. Paine and Kittel (1958) and Loomis and Gerber (1968) both investigated dimensional change of single-crystal uranium under irradiation. Their findings were in general agreement with one another, and showed a contraction of the bulk crystal in the a direction, an expansion in the b direction, and no change in the c direction. While the AIMD results of point defect-induced strains presented here agree with the lack of dimensional change in the c direction, exactly the opposite behavior is observed in Figure 7 for the a and b directions. This work showed that interstitials are the dominant strain inducers, producing expansion in the a direction and contraction in the b direction. Due to the relatively large induced strains from interstitials, it is expected that these types of point defects would display strongly anisotropic diffusion. Given that isolated point defects lead to the opposite trends from those observed experimentally, it is possible that defect clusters or defect loops display different strain behaviors than individual point defects, and are the features responsible for irradiation growth in α-U. Such research is beyond the scope of AIMD, which is currently sufficiently computationally expensive to exclude the possibility of analysis of extended defects. Thus, classical molecular dynamics simulations can be performed and validated against the work presented here. Ideally, in addition to further computational studies, the experimental results from over 60 years ago should be reproduced and modern characterization methods applied to extract fundamental defect populations, densities, and orientations.

In this study, AIMD simulations were performed to calculate temperature-dependent properties of the α phase of U from 100 to 800 K. The thermal expansion of each lattice constant was analyzed and compared to the experimental literature. Although the general direction of thermal expansion corresponded to the experimental findings (positive in a and c, negative in b), the magnitudes of coefficient of thermal expansion for the a and b directions was greater than that observed experimentally. However, the volumetric expansion corresponded nearly exactly to the experimental findings in the literature. The heat capacity was calculated as a function of temperature and showed a general agreement with experimental results, although the computationally determined heat capacity was slightly less sensitive to temperature. Variation in point defect formation energies was determined, with an average interstitial formation energy of 3.40 eV, and an average vacancy formation energy of 1.46 eV, over the entire temperature range. Lattice strains resulting from point defects were determined, showing that interstitials induce significantly more strain than vacancies, and the nature of that strain is highly dependent on the individual lattice directions, where the a direction expands and the b contracts due to the introduction of an interstitial. Finally, all of these results seem to indicate a two-region behavior, with the transition at 400 K. The thermal expansion/contraction in the a/b direction occurs rapidly from 100 up to 400 K, after which a more linear and gradual change in the lattice constant takes place. Defect formation energies decrease in magnitude up to 400 K, above which they tend to slightly increase. Defect lattice strains for interstitials decrease/increase in the a/b direction up to 400 K, above which they maintain essentially the same magnitude. Coordination and electronic structure investigations provided no insight into the possible cause of this potential transition in behavior.

This work has shown that isolated point defects cannot be the primary driving force responsible for the significant irradiation-induced growth of α-U observed experimentally. This work can be utilized as the basis for mesoscale simulations that take into account defect populations, defect diffusion, and defect strain behaviors. It is emphasized that additional computational work utilizing classical molecular dynamics on extended defects should be performed, in addition to experimental irradiation and characterization, to further elucidate the underlying phenomena governing the interactions of point defects with the α-U crystal structure.

Figures, tables, and images will be published under a Creative Commons CC-BY license and permission must be obtained for use of copyrighted material from other sources (including re-published/adapted/modified/partial figures and images from the internet). It is the responsibility of the authors to acquire the licenses, to follow any citation instructions requested by third-party rights holders, and cover any supplementary charges.