Mechanical and Thermal Properties for Uranium and U–6Nb Alloy from First-Principles Theory

: Elasticity, lattice dynamics, and thermal expansion for uranium and U–6Nb alloy (elastic moduli) are calculated from density functional theory that is extended to include orbital polarization (DFT+OP). Introducing 12.5 at.% of niobium, substitutionally, in uranium softens all the c ii elastic moduli, resulting in a signiﬁcantly softer shear modulus (G). Combined with a nearly invariant bulk modulus (B), the quotient B/G increases dramatically for U–6Nb, suggesting a more ductile material. Lattice dynamics from a harmonic model coupled with a DFT+OP electronic structure is applied for α uranium, and the obtained phonon density of states compares well with inelastic neutron-scattering measurements. The Debye temperature associated with the lattice dynamics falls within the range of experimentally observed Debye temperatures and it also validates our quasi-harmonic (QH) phonon model. The QH Debye–Grüneisen phonon method is combined with a DFT+OP electronic structure and used to explore the anisotropic thermal expansion in α uranium. The anomalous negative thermal expansion (contraction) of the b lattice parameter of the α -phase orthorhombic cell is relatively well reproduced from a free-energy model consisting of QH-phonon and DFT+OP electronic structure contributions.


Introduction
Uranium is an essential metal that has many different applications. Of course, its nuclear attributes are important, but its material properties are also very relevant. It is a heavy (high-density) metal that is mechanically strong relative to other nuclear materials but is prone to oxidation (corrosion). Uranium's mechanical properties are to a degree a consequence of its rather unusual crystal structure (orthorhombic) that induces greater anisotropy than materials with higher symmetry structures, such as the d-transition metals. The material properties and the anisotropy are reflected in the elastic moduli that represents the character and strength of the atomic bonding, that in the case of uranium is dominated by delocalized and bonding 5f electrons [1]. The anisotropy is also revealed in an anomalous and poorly understood negative thermal expansion (contraction) of one of its orthorhombic lattice constants.
Adding solid-solution elements from the series of d-transition metals has a beneficial effect on corrosion resistance and improves the mechanical properties. Typically, these alloying elements are soluble at higher temperatures in the body-centered cubic (bcc) phase, γ-uranium, but they can be retained in other phases through quenching [2], as we will discuss below. It has been shown that for U-Nb-specific alloys, the alloying improves the mechanical properties, including ductility, relative to pure uranium [3,4]. Additionally, the addition of niobium to uranium is particularly beneficial, as it significantly increases oxidation resistance because the niobium metal is distributed uniformly in the uranium matrix [5]. These uranium-niobium alloys furthermore exhibit a shape memory effect (SME) that depends on temperature and strain conditions.

Computational Methods
All results presented in this article are derived from the calculated free energy of uranium or niobium atoms and their geometrical arrangement. The theoretical framework is established by DFT, which only fundamental approximation is that of the electron exchange and correlation. There are many different alternatives for this approximation that govern the quantum-mechanical interaction between the electrons. They make different assumptions that are better or worse depending on the material in question. Fortunately, the generalized gradient approximation (GGA) [34] is known [35] to be quite good for actinide metals, including uranium. There are better variants for niobium (a 4d transition metal) but the niobium content in the studied uranium-niobium alloy is relatively low (no more than 12.5% of the atoms) and therefore the GGA remains a good choice.
Because the GGA depends on the spin-polarized electron density, the effects of magnetism and the relativistic spin-orbit interaction are included in the theory. However, in the conventional DFT formulation, there are no electron currents, so orbital-orbital interaction is ignored. Nonetheless, for the actinides and particularly plutonium [36], this interaction is important and there are approximate ways to amend DFT to include it. The "orbital polarization" (OP) scheme defines a parameter-free self-consistent procedure that is computationally tractable [37] and we employ it here (DFT+OP), even though its effect in uranium is less pronounced than in plutonium [38].
For best accuracy, it is essential to minimize numerical or other approximations beyond the GGA, and the all-electron technique is preferred over methods that rely on pseudopotential or other approximations for core electrons not included in the valence. Here, we utilize the all-electron full-potential linear muffin-tin orbitals (FPLMTO) method [39]. Basis functions, electron densities, and potentials are calculated with no geometrical approximation, and these are expanded in spherical harmonics (with a cutoff l max = 8) inside non-overlapping (muffin-tin) spheres around the atoms and in Fourier series in the space outside these muffin-tin spheres. The muffin-tin sphere radius is fixed and close to 0.8 of that of a sphere with a volume equal to the atomic volume. All calculations utilize semi-core states and valence states with two fixed energy parameters each for the s semi-core state, p semi-core state, and the valence states. The number of k points included in the electronic structure calculations depends on the crystal structure, but we generally use~2000 k points or more for one atom/cell calculations, and less for structures with more atoms. Each energy eigenvalue is broadened with a Fermi-Dirac function at room temperature (300 K) unless free energies at higher temperatures are considered.
We introduce a 16-atom supercell approximation for the U-6Nb alloy by replacing 2 of the 16 uranium atoms with niobium in an ordered U 14 Nb 2 compound, either in the monoclinic α" phase or the orthorhombic α phase. This assumption is analogous to the study by Wang et al. [32]. Geometrical energy optimization (structural relaxation) is performed for the genuine phases (α for uranium and α" for the U-6Nb alloy). For α uranium, the relaxation produces axial ratios b/a = 2.045, c/a = 1.755, an internal parameter, y = 0.1009, and an atomic volume 20.67 Å 3 . Correspondingly, the uranium-niobium alloy in the α" phase has axial ratios b/a = 2.022, c/a = 1.753, monoclinic angle 91.3 • , internal parameters y 1 = 5/6 and y 2 = 1/3, and atomic volume 20.48 Å 3 .
The thirteen (nine) elastic constants for the monoclinic (orthorhombic) α" (α) structure are calculated for the pure element as well as the uranium-niobium alloy. The procedure is to apply about 20 strains (within ±1%) [40,41] corresponding to the various moduli, leastsquare fit the free-energy response to a 4th order polynomial and extract the coefficient for the 2nd order term [40,41]. No structural relaxations are considered during these small strains. To calculate the shear (G V ) and bulk (B V ) moduli, we use standard Voigt expressions. The ratio B V /G V indicates a materials brittleness (small value) or ductility (large value) [42].
The phonon density of states (DOS) is calculated in the harmonic approximation for α uranium. We employ the small-displacement method, utilizing a 36-atom supercell and two small atomic displacements, as implemented in the PHONOPY code [43]. To improve on the accuracy of the phonons, we performed an additional supercell calculation without displacements so that we could subtract any residuals in the force-matrix result. Technically, the forces are determined numerically from an approach that relies on moving each atom and calculating the energy response. This somewhat cumbersome procedure is robust and generates reliable forces from an electronic structure that is significantly perturbed by spin-orbit coupling and orbital polarization. Under these circumstances, the atomic forces are not easily computed accurately using standard methods; see previous discussions [44]. The Debye temperature is obtained from integrating over the computed phonon density of states in the conventional fashion.
The thermal expansion for each of the a, b, and c lattice constants of α uranium is derived from varying the b/a axial ratio and evaluating the free energy from 300 K to 900 K. It is known [45,46] that the a and c lattice parameters have nearly identical thermal expansions above 300 K, thus motivating us to simplify the problem and keep c/a constant (1.755).
In our approach, the free energy consists of an electronic contribution, obtained directly from our electronic structure, and a lattice-vibration term. The electron-phonon coupling is not easily formulated, and because of the relatively low temperatures, this interaction is assumed to be negligible. The phonon contribution is expressed within Debye-Grüneisen (DG) quasi-harmonic (QH) theory and computed utilizing our own codes [47]. The DG approach is validated by comparing the Debye temperature (θ D ) with that obtained independently from our first-principles phonon density of states (both θ D = 240 K). From the calculated free-energy minimum we determine the b/a axial ratio as a function of temperature. The absolute thermal expansions for the a, b, and c lattice constants are then established from the total-volume expansion acquired from the QH-DG model.

Elastic Constants
Elastic constants reflect important material properties related to the strength and character of the interatomic bonding. They directly connect to macroscopic measurables such as ductility, brittleness, and resistance to shear and pressure. Furthermore, elasticity is an essential ingredient in constitutive modeling and involved in the shape memory effect that exists in the U-Nb materials. Elastic moduli can also be used to constrain interatomic potentials that are developed and employed in molecular dynamics simulations.
Here, we focus on the uranium-niobium alloy with about 6 wt.% niobium (~14 at.% Nb), often denoted U-6Nb (approximated here by an ordered U 14 Nb 2 compound). Because we are interested in the influence of niobium on the elastic properties, we explore unalloyed uranium as well, and include its ground-state orthorhombic α phase in the study as a reference in addition to the monoclinic α" (U-6Nb) phase.
The calculated elastic constants for uranium and U 14 Nb 2 in the α" and α phases are listed in Table 1. We are not aware of either calculations or measurements for the uraniumniobium alloy, but for α uranium the data are consistent with others' modeling [40,[48][49][50][51] and rather close to the experimental data [52], particularly when extrapolated to zero temperature. The good comparison between experiment and modeling for α uranium gives us confidence that the theoretical framework is also sound and appropriate for the uranium-niobium alloy. Table 1. The independent elastic moduli (Mbar) for the α" and α phases of uranium and uraniumniobium alloy. G V and B V are the corresponding shear-and bulk-modulus Voigt averages. The quotient B V /G V (Pugh ratio) relates to malleability or ductility (high value) or brittleness (low value) [42]. For α-U, experimental data [52] shown in parentheses are extrapolations to zero temperature. Note, α" and α are hypothetical phases for uranium and uranium-niobium alloy, respectively, shown here for comparison. The obtained elastic moduli are not dramatically different between U 14 Nb 2 and uranium in either the α or α" phase, and that is likely because only 12.5% of the uranium atoms are replaced by niobium atoms. However, we find a general trend of substantial elastic softening in both the α and α" phase. All c ii moduli are significantly smaller for the alloy containing niobium. The bulk modulus, however, is relatively uniform between uranium and the uranium-niobium alloy in both phases. The shear modulus (G V ), on the other hand, is notably smaller for the alloy. Consequently, the Pugh ratio B V /G V increases by an appreciable 30-40% for the alloy, suggesting a great improvement in ductility over the pure element.
We also performed reference calculations with 1/16 niobium (U 15 Nb) in the α phase only, and in this case the ductility (in terms of the Pugh ratio) increased by about 17% relative to elemental uranium. There is a proportionality between niobium content and ductility, at least for smaller amounts of niobium. The reason for the better ductility is directly traced to the softening of the c ii single-crystal elastic constants that we compute for the alloy (Table 1).
Another interesting observation in Table 1 is that c 11 , c 22 , and c 33 are quite different for the studied phases. These three moduli correspond to the energy cost to independently strain along the a, b, and c axis, respectively. Consequently, there is a great deal of anisotropy in the crystal both in the α and α" phases for uranium and the uraniumniobium alloy. Particularly, c 22 stands out and is smaller than the other two. This is correlated with the anisotropic behavior of the b lattice constant with temperature, and we shall return to this in Section 3.3.
Finally, the present U 14 Nb 2 elastic constants in the α" and α phases have practical value because they can help constrain interatomic potentials or guide constitutive modeling for the shape memory U-6Nb alloy.

Lattice Dynamics
In addition to the elasticity, lattice vibrations or phonons provide important thermal and structural information about the material. For example, soft phonons could signal structural phase transitions, while stiffer phonons usually imply greater elastic constants and strength. The phonons also contain temperature-dependent information and parameters for thermodynamic modeling, such as the Debye temperature or the Grüneisen parameter.
By relating the theoretical phonons with experimental observations, one can assess the accuracy and reliability of the model. If a material's calculated atomic density, crystal structure, elastic constants, and phonons all compare favorably with measurements, then the chemical bonding predicted by the model must be rather accurate as well. Recently, excellent phonons were calculated for a far more controversial material than α uranium, namely α plutonium [53], thus validating the electronic structure model assumed for plutonium [36].
Our approach for modeling the α-uranium phonon density of states is the same as that recently performed for α plutonium [53]. One difference is that for plutonium, non-negligible spin and orbital magnetic moments are formed [54], while for α uranium, the electrons do not spontaneously spin polarize. Nonetheless, effects of both spin-orbit interaction and orbital polarization are present and accounted for in our theory for uranium.
In Figure 1, we show our DFT+OP phonon density of states together with experimental data [55] and least-square fits of force-constant models to experimental data [56]. In this figure, the notations are repeated from [55], and for more details we refer to that publication. The neutron-scattering data [55] show significant scatter, as do the force-constant models, but present in all data is the two-peak feature with one peak at 10 meV and the other around 12.5-14.5 meV. The DFT+OP model accurately reproduces this double-peak feature, while at the highest energies there is a shift between DFT+OP on one hand and the neutron data on the other, with the force-constant models in the middle.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 6 of 12 and strength. The phonons also contain temperature-dependent information and parameters for thermodynamic modeling, such as the Debye temperature or the Grüneisen parameter. By relating the theoretical phonons with experimental observations, one can assess the accuracy and reliability of the model. If a material's calculated atomic density, crystal structure, elastic constants, and phonons all compare favorably with measurements, then the chemical bonding predicted by the model must be rather accurate as well. Recently, excellent phonons were calculated for a far more controversial material than α uranium, namely α plutonium [53], thus validating the electronic structure model assumed for plutonium [36].
Our approach for modeling the α-uranium phonon density of states is the same as that recently performed for α plutonium [53]. One difference is that for plutonium, nonnegligible spin and orbital magnetic moments are formed [54], while for α uranium, the electrons do not spontaneously spin polarize. Nonetheless, effects of both spin-orbit interaction and orbital polarization are present and accounted for in our theory for uranium.
In Figure 1, we show our DFT+OP phonon density of states together with experimental data [55] and least-square fits of force-constant models to experimental data [56]. In this figure, the notations are repeated from [55], and for more details we refer to that publication. The neutron-scattering data [55] show significant scatter, as do the force-constant models, but present in all data is the two-peak feature with one peak at 10 meV and the other around 12.5-14.5 meV. The DFT+OP model accurately reproduces this doublepeak feature, while at the highest energies there is a shift between DFT+OP on one hand and the neutron data on the other, with the force-constant models in the middle. The discrepancy at the highest energies may be due to the lack of explicit temperature dependence of the phonons in the method, or that only harmonic lattice vibrations were considered in the model. One must also acknowledge that the experimental data may have rather large error bars, particularly at the highest energies (>12 meV) where the data points are more scattered; see Figure 1.
The Debye temperature (θD) can be expressed in terms of a Brillouin zone average of the phonon frequencies [57], and from the DFT+OP phonon DOS we determine θD = 240 K. Experimentally, there is a wide range (183-259 K) of room-temperature Debye temperatures reported and listed in [55] and [58], with the 11 data-point average being θD = 217 The discrepancy at the highest energies may be due to the lack of explicit temperature dependence of the phonons in the method, or that only harmonic lattice vibrations were considered in the model. One must also acknowledge that the experimental data may have rather large error bars, particularly at the highest energies (>12 meV) where the data points are more scattered; see Figure 1.
The Debye temperature (θ D ) can be expressed in terms of a Brillouin zone average of the phonon frequencies [57], and from the DFT+OP phonon DOS we determine θ D = 240 K. Experimentally, there is a wide range (183-259 K) of room-temperature Debye temperatures reported and listed in [55] and [58], with the 11 data-point average being θ D = 217 K. The Debye temperature is also known to be rather temperature dependent and decreases significantly with temperature [55].
We mention that another theoretical (DFT) study of α uranium produced lattice dynamics in good agreement with neutron-scattering data as well [48]. In that study, the focus was on the low-temperature properties of α uranium and specifically the distortions associated with the charge-density wave.

Thermal Expansion
The thermal expansion of α uranium is strongly anisotropic [45,46]. Above room temperature, two (a and c) of the three lattice constants display a normal (and quite similar) positive expansion with temperature, while b anomalously shrinks or contracts with a negative thermal expansion. The volumetric thermal expansion, however, remains positive above room temperature. Negative thermal expansion is unusual in general, and among the actinides the δ-phase of plutonium is probably best known for this exceptional behavior. For δ plutonium, its origin has been proposed to be thermally excited longitudinal spin fluctuations [36]. For α uranium, the formation of magnetic moments has been suggested [59] as an explanation for the anisotropic behavior. To our knowledge, a fundamental and detailed understanding of this phenomenon in α uranium has not yet been provided. Regardless of cause, the expansion behavior of uranium is important since it relates to the thermal dependency of the mechanical properties and the shape memory effect in the uranium-niobium alloy.
Because magnetic moments do not survive (collapse) in α uranium in our firstprinciples model, any negative thermal expansion needs to have a different explanation than magnetism for the model to be correct. To explore a more straightforward origin, namely, a combination of conventional electronic and lattice-vibration thermal effects on the crystal structure (axial ratios), we study the free energy as a function of the b/a axial ratio and temperature. We simplify the model by keeping the c/a axial ratio constant at 1.755 and the internal parameter y at 0.1009. This is a reasonable simplification because the y parameter is rather insensitive to small changes in the lattice constants (not shown) and the a and c lattice constant expansions are nearly identical [45,46] not too far above room temperature.
Quasi-harmonic theory in the framework of the Debye-Grüneisen model is employed for the lattice-vibrational contribution to the free energy, and the electronic contribution is computed within the DFT+OP approach. The calculated quasi-harmonic Debye temperature is identical (240 K) to that obtained independently from our DFT+OP phonon calculation, and that reassures us of the appropriateness of the Debye-Grüneisen model for the vibrational contribution.
In Figure 2, we show free energies, relative to the ground-state α-uranium phase (b/a = 2.045, c/a = 1.755, and y = 0.1009), for structures with different b/a axial ratios as functions of temperature. The top and middle panel display the lattice and electronic freeenergy contributions, respectively. We note that increasing temperature favors a decrease in the b/a axial ratio for both the lattice and electronic contributions to the free energy. These quantities are similar in magnitude and equally important for the net result. The electronic term at T = 0 K reflects the structural energetics at zero temperature and shows a weak energy dependence on the b/a axial ratio (~0.25 mRy/atom between b/a = 2.045 and b/a = 1.995). This weak energy dependence is essential because otherwise the b/a ratio would be insensitive to an increase in temperature (and it obviously is not). The bottom panel of Figure 2 shows the total free energy (Ftotal = Felectronic + Flattice) and one finds that the lowest energies are predicted to occur for different b/a ratios at different temperatures. First, we have b/a = 2.045 at zero temperature, and then we see transitions to (a) 2.035, (b) 2.025, (c) 2.015, (d) 2.005, and finally (e) 1.995 values of b/a with increasing temperatures. These transitions occur only because of the delicate balance between the structural energetics at zero temperature and the temperature-dependent lattice and electronic free-energy terms. The net effect of these contributions is a decrease in the b/a axial ratio with temperature for α uranium. Because we know the total volumetric thermal expansion from the Debye-Grüneisen model, that incidentally agrees rather well with the experiment [47], we can deduce the separate a, b, and c lattice constants' temperature dependence.
In Figure 3, we plot the a, b, and c lattice constants as functions of temperature together with experimental data [45]. At room temperature, these lattice parameters are about 1% or less in error relative to the experiment, and the thermal dependence for all three, including the anomalous behavior for the b parameter, is also reproduced quite well. The bottom panel of Figure 2 shows the total free energy (F total = F electronic + F lattice ) and one finds that the lowest energies are predicted to occur for different b/a ratios at different temperatures. First, we have b/a = 2.045 at zero temperature, and then we see transitions to (a) 2.035, (b) 2.025, (c) 2.015, (d) 2.005, and finally (e) 1.995 values of b/a with increasing temperatures. These transitions occur only because of the delicate balance between the structural energetics at zero temperature and the temperature-dependent lattice and electronic free-energy terms. The net effect of these contributions is a decrease in the b/a axial ratio with temperature for α uranium. Because we know the total volumetric thermal expansion from the Debye-Grüneisen model, that incidentally agrees rather well with the experiment [47], we can deduce the separate a, b, and c lattice constants' temperature dependence.
In Figure 3, we plot the a, b, and c lattice constants as functions of temperature together with experimental data [45]. At room temperature, these lattice parameters are about 1% or less in error relative to the experiment, and the thermal dependence for all three, including the anomalous behavior for the b parameter, is also reproduced quite well. While the trend is good, we notice there is some scatter in the low-temperature results, and we attribute that to the constraint of using discrete axial b/a ratios. Only at the highest temperatures is the predicted trend exaggerated. This may not be surprising because the quasi-harmonic model is expected to be less accurate at higher temperatures, where anharmonic lattice vibrations may become non-negligible. In addition, upon closer inspection of the experiments, we see that the c lattice constant has a somewhat greater thermal expansion than a at high temperatures. This behavior is ignored in the model because we assume a frozen c/a axial ratio. In a more advanced treatment, where c/a is allowed to change with temperature, we expect an even better agreement with the experiment for the b/a ratio, because at high temperatures, the a lattice constant is slightly underestimated, resulting in an excessive decrease in the b lattice constant.

Discussion
We have studied uranium metal and U-6Nb alloy from the first-principles theory. For the uranium-niobium alloy, we focus on elastic moduli, and the essential information on its elasticity has not been reported until now. The elastic properties are important in terms of modeling, i.e., they can constrain interatomic potentials for molecular dynamics simulations, or they can explicitly be used in constitutive modeling. The elastic constants also explain the appreciated fact that the uranium-niobium alloys are more ductile than the pure metal because of a significant alloy-driven softening of the cii moduli. For both the metal and the alloy (α or α" phase), we find that c11, c22, and c33 are very different. In particular, c22, that correlates with a change in the b lattice parameter, is significantly smaller than the other two. This crystal anisotropy shows up dramatically in the thermalexpansion behavior as well. While the trend is good, we notice there is some scatter in the low-temperature results, and we attribute that to the constraint of using discrete axial b/a ratios. Only at the highest temperatures is the predicted trend exaggerated. This may not be surprising because the quasi-harmonic model is expected to be less accurate at higher temperatures, where anharmonic lattice vibrations may become non-negligible. In addition, upon closer inspection of the experiments, we see that the c lattice constant has a somewhat greater thermal expansion than a at high temperatures. This behavior is ignored in the model because we assume a frozen c/a axial ratio. In a more advanced treatment, where c/a is allowed to change with temperature, we expect an even better agreement with the experiment for the b/a ratio, because at high temperatures, the a lattice constant is slightly underestimated, resulting in an excessive decrease in the b lattice constant.

Discussion
We have studied uranium metal and U-6Nb alloy from the first-principles theory. For the uranium-niobium alloy, we focus on elastic moduli, and the essential information on its elasticity has not been reported until now. The elastic properties are important in terms of modeling, i.e., they can constrain interatomic potentials for molecular dynamics simulations, or they can explicitly be used in constitutive modeling. The elastic constants also explain the appreciated fact that the uranium-niobium alloys are more ductile than the pure metal because of a significant alloy-driven softening of the c ii moduli. For both the metal and the alloy (α or α" phase), we find that c 11 , c 22 , and c 33 are very different. In particular, c 22 , that correlates with a change in the b lattice parameter, is significantly smaller than the other two. This crystal anisotropy shows up dramatically in the thermal-expansion behavior as well.
The phonon density of states is calculated for α uranium and compares well with neutron scattering. This fact, and that theoretical elastic moduli are in good agreement with measured data (particularly when compared at the same 0 K temperature), gives us confidence in the DFT+OP model used here for uranium and the uranium-niobium alloy. The Debye temperature obtained from the DFT+OP phonon DOS is consistent with that of the quasi-harmonic model, thus validating the latter as a sound approach for the thermal-expansion study.
For α uranium, we investigate the unusual thermal-expansion behavior. We show that the anisotropic and negative thermal expansion does not need to originate from exotic (magnetic) electronic structure. Instead, a delicate balance between zero-temperature energetic and temperature-dependent electronic and lattice-vibration contributions to the free energy fully explains the anomaly. In this approach, the electronic term is determined from DFT+OP electronic structure with a Fermi-Dirac broadening associated with the electronic density of states, and the phonon term is based on the Debye-Grüneisen quasi-harmonic scheme. The agreement between the observed anomalous lattice-constant anisotropy and that predicted by the model is satisfactory, at least below the highest temperatures.