Assessing Relativistic E ﬀ ects and Electron Correlation in the Actinide Metals Th to Pu

: Density functional theory (DFT) calculations are employed to explore and assess the e ﬀ ects of the relativistic spin–orbit interaction and electron correlations in the actinide elements. Speciﬁcally, we address electron correlations in terms of an intra-atomic Coulomb interaction with a Hubbard U parameter (DFT + U ). Contrary to recent beliefs, we show that for the ground-state properties of the light actinide elements Th to Pu, the DFT + U makes its best predictions for U = 0. Actually, our modeling suggests that the most popular DFT + U formulation leads to the wrong ground-state phase for plutonium. Instead, extending DFT and the generalized gradient approximation (GGA) with orbital–orbital interaction (orbital polarization; OP) is the most accurate approach. We believe the confusion in the literature on the subject mostly originates from incorrectly accounting for the spin–orbit (SO) interaction for the p 1 / 2 state, which is not treated in any of the widely used pseudopotential plane-wave codes. Here, we show that for the actinides it su ﬃ ces to simply discard the SO coupling for the p states for excellent accuracy. We thus describe a formalism within the projector-augmented-wave (PAW) scheme that allows for spin–orbit coupling, orbital polarization, and non-collinear magnetism, while retaining an e ﬃ cient calculation of Hellmann–Feynman forces. We present results of the ground-state phases of all the light actinide metals (Th to Pu). Furthermore, we conclude that the contribution from OP is generally small, but substantial in plutonium.


Introduction
In the last several decades, various computational techniques to calculate the electronic structure and ground-state properties of actinides have been suggested in the literature . In these studies, a comparison of the calculated equilibrium volumes with experiments has often served the purpose to evaluate the quality of the modeling and the underlying physics of the 5f electrons. From these works a consistent picture emerges for thorium through neptunium where the 5f electrons can be described as bonding and delocalized (band) electrons. The implication is that the actinide volumes decrease in a parabolic fashion with atomic number for Th to Pu, reminiscent of the non-magnetic d-transition metals. The physical interpretation, given by Friedel [27], is that for the first half of the series, bonding electrons contract the volumes while for the second half, anti-bonding electrons expand the lattice. The parabolic behavior is readily verified by assuming a rectangular electron density of states and the corresponding electron-energy dependence on the volume [27].
Plutonium, on the other hand, deviates slightly from this simple trend. In this case, relativistic effects [4], crystal structure [8,10], and possible electron correlation effects beyond the local density approximation [9] are suggested as explanations for some of the discrepancy noted between theory and experiments. However, a much-improved description of plutonium is obtained, still within the 5f-band picture, when magnetic effects are addressed [13]. Acknowledging these magnetic interactions

Computational Methods
We show results derived from density functional theory that was developed in the mid-1960s [42,43]. Density functional theory is actually correct, but it relies on practical approximations. The central assumption relates to the electrons exchange and correlation interactions. For the actinide metals, we employ the generalized gradient approximation that is more proper than the previously suggested local density approximation [6]. Specifically, we are applying the Perdew-Burke-Ernzerhof (PBE) parameterization of the generalized gradient approximation (GGA) [44]. To our knowledge, there is no established, parameter-free, functional that results in larger atomic volumes in the actinides than that of the GGA. Therefore, there is no choice to replace the GGA to expand the volumes, while a +U approach can accomplish this objective and has often been applied for the actinides, as discussed in the introduction.
The actinide metals exhibit interesting crystal structures, beginning with an innocent face-centered cubic for Th with a gradual escalation of complexity to tetragonal (Pa), orthorhombic (U and Np), and finally culminating with a monoclinic structure (α-plutonium). These structures are described in the compilation by Donohue [45]. For best possible accuracy, when calculating the equilibrium volumes, we have optimized these crystal structures (not the cubic phase) so that the model gives the lowest total energy (structural relaxation). This relaxation has generally a rather weak effect on the actinide volumes, except for α-plutonium where it is significant. In this instance, a structural relaxation contracts the volume substantially relative to a calculation where the experimental (room-temperature) structure is assumed [40,46].
It has been reported that orbital-orbital interactions, absent in conventional DFT, is important for plutonium [14,40,47]. This interaction is efficiently included in an orbital-polarization scheme for some of the calculations presented in this study. In terms of the technical details of the computations, we are applying three techniques. The bulk of the calculations are performed with the pseudopotential plane-wave Vienna ab initio simulation package (VASP) code.
In VASP, the computational scheme relies on the GGA with PBE exchange and correlation functional [44] and the projector-augmented-wave (PAW) method with plane-wave basis set as implemented in VASP [48][49][50]. The computational parameters consist of an energy cut-off of 450 eV, energy convergence of 100 eV, and a Monkhorst-Pack scheme k-point mesh of typically 12 × 12 × 12 for a 1-atom cell and correspondingly less for multi-atom cells. In the cases where we perform DFT + U computations, we employ the technique suggested by Dudarev et al. [51] for the 5f electrons with an effective U = 0.4 eV.
The spin-orbit interaction and orbital polarization are included, allowing for non-collinear spins. For the spin-orbit coupling, we have introduced a correction that retains the scalar-relativistic treatment for the p states. Orbital polarization has been implemented similarly to the scheme proposed by Eriksson et al. [52]. In this approach, the DFT total-energy functional, E DFT , is amended by a quadratic term favoring orbital polarization, where E 3 I is the Racah parameter that is a linear combination of the radial Slater integrals F 2 , F 4 , and F 6 over f-electron orbitals within a muffin-tin sphere centered on site I. For consistency and simplicity, we use a constant Racah parameter that is an average of spin-decomposed values obtained from all-electron calculations (see below). L f I,z is the z component of the angular-momentum operator within the same sphere operating on the f states (5f in this case). This formulation of the orbital-polarization Hamiltonian has mainly been implemented in conjunction with spin-orbit coupling in all-electron codes within the second-variation approximation assuming collinear spin polarization.
Analogous to the conventional formulation of the orbital polarization (OP) method, while having the computational expediency of plane-wave codes and accurate Hellman-Feynman forces, we make use of a variational total-energy formulation of the OP Hamiltonian in the context of non-collinear magnetism. It is easily constructed within the PAW method following the DFT + U implementation by Bengone et al. [53]. See Appendix A for a brief technical description of the PAW method for the purpose of implementation of the orbital-polarization Hamiltonian.
As mentioned, Equation (1) is proposed for collinear magnetism with the spin moment pointing in the z direction. This formalism can be generalized to non-collinear magnetism with arbitrary spin orientations by modifying the quadratic energy functional within each muffin-tin radius to favor orbital polarization in the direction of the expectation value of the local self-consistent spin vector. It is then straightforward to propose the following total-energy functional, consistent with a non-collinear generalization of Equation (1), where L  In addition to our plane-wave calculations, we apply two all-electron methods to compare our pseudopotential results with. First, a full-potential linear muffin-tin orbital (FPLMTO) method [54] that has been detailed [55]. Here, no approximations are applied for the core electrons, unlike the pseudopotential method. The core-electron approximation made in the latter technique is computationally expedient but will cause inaccuracies that are not present in all-electron calculations.
The present implementation does not make any assumptions beyond the GGA. Basis functions, electron densities, and potentials are calculated without any geometrical approximation, and these are expanded in spherical harmonics (with a cut-off L max = 8) inside non-overlapping (muffin-tin, MT) spheres surrounding each atom and in Fourier series in the region between these muffin-tin spheres. One has to define an MT sphere with a radius, s MT , and here it is chosen so that s MT /s WS~0 .8, where s WS is the Wigner-Seitz (atomic sphere) radius. The radial parts of the basis functions inside the MT spheres are calculated from a wave equation for the L = 0 component of the potential that includes all relativistic corrections including spin-orbit coupling for d and f states, but not for the p states [20]. It should be noted that excluding the p states also removes the strong and artificial dependence of the calculated volumes on the s MT /s WS ratio [2] that is also not present in the exact relativistic treatment (see below).
The orbital polarization is only considered for d and f states. In the actinides, however, the OP contribution from the d states is minimal and ignored in the VASP calculations. The formulation is self-consistent and parameter-free and given in [52]. The self-consistent FPLMTO Racah parameters for the 5f states are employed in the VASP calculations described above. All crystal-structure geometries were relaxed and for neptunium and plutonium, that have more complex phases, we apply a scheme [56] that avoids the necessity of Hellman-Feynman forces that are troublesome to compute.
The number of k points included in the calculations depends on the particular crystal structure, but we generally use~1000 k points or more for one atom/cell calculations and fewer for cells with many atoms. Each energy eigenvalue is broadened with a Gaussian having a width of 20 mRy.
Second, we employ a full-potential linear augmented plane-wave (FPLAPW) method as implemented in FlapwMBPT (many body perturbation theory) [57]. This approach has many similarities with the FPLMTO method introduced above. As in the FPLMTO method, the space is divided into non-overlapping MT spheres and an interstitial region. The electronic density and effective potential are expanded in spherical harmonics in the MT spheres and in plane waves in the interstitial region (exactly as in the FPLMTO method).
One conceptional difference between FLAPW and the other two methods is that orbital polarization is not applied in the FLAPW calculations. For this reason, Th, Pa, and U were studied with FLAPW, but not Np and Pu where orbital polarization is more important. Another conceptional difference is that FlapwMBPT allows for a complete inclusion of relativistic effects through the four-spinor Dirac equation. Particularly important is that the spin-orbit interaction is here formulated exactly.
Other technical differences between FLAPW and FPLMTO are mostly related to the use of different basis functions in the interstitial region (plane-waves and Hankel functions, respectively) and in the augmentation of them with the functions inside the MT spheres. Details of the implementation of the FLAPW method can be found in [58]. Specific values of the different cut-offs used in the FLAPW calculations are summarized in Table 1. The crystal-structure parameters given in Table 1 are obtained from structural relaxations using the FPLMTO method.

Results
The formally correct way to include spin-orbit interaction in an electronic-structure calculation is to solve the four-spinor Dirac equation. Here, we have applied a FPLAPW method that indeed does this (FlapwMBPT [57]), and we regard the corresponding results as the gold standard for our comparisons. We mentioned [20] that utilizing a scalar-relativistic basis set causes inaccuracies for the p states, but it has been suggested that this problem can be circumvented for the actinides by ignoring spin-orbit interaction for these states [20]. This realization has not been widely recognized in studies using plane-wave pseudopotential techniques, and we therefore reiterate this fact by showing our results in Table 2.  Table 2 presents the effect of spin-orbit coupling on the volumes of the actinide elements. For Th, Pa, and U, FPLAW-Dirac calculations provide the reference, while for Np and Pu, that have substantially more complex and large crystal structures, we show results from VASP and FPLMTO. For VASP and FPLMTO the perturbative form of SO is included but for the p states. All three methods give essentially identical results for Th, Pa, and U, suggesting that ignoring spin-orbit interaction for the p states is a very good approximation. The same conclusion was reached in previous studies [59,60] where comparisons were made with the exact muffin-tin orbital method that also solves the four-spinor Dirac equation.
The trend is that spin-orbit coupling is of increasing importance proceeding from thorium to neptunium, causing an expansion of the volumes ( Table 2). These metals are all non-magnetic and the explanation for this behavior is straightforward. As mentioned above, the atomic volumes for Th to Np decrease with atomic number in roughly a parabolic fashion that is due to a gradual population of bonding 5f states. No anti-bonding 5f electrons exist because Th to Np all have less than half-filled 5f band. Inclusion of spin-orbit interaction, however, causes some anti-bonding 5f states to be occupied, and these tend to expand the atomic volumes. The reason is that anti-bonding states mix with bonding states as the 5f band is divided into overlapping 5f 5/2 and 5f 7/2 sub-bands that comprise both bonding and anti-bonding states. The result of this is a reduction of the effective number of bonding electrons and an increase in the atomic volumes (Table 2).
For plutonium, on the other hand, the situation becomes more complicated because here both spin and orbital polarization play important roles as well. In addition, for α-plutonium the crystal and electronic structure correlate sensitively so that the atomic volume cannot be obtained simply from extrapolating the trend established by Th to Np. For plutonium, one must thus carry out accurate electronic-structure calculations, recognizing the magnetic interactions and relaxing the crystal structure, to predict the atomic volume accurately.
Next, we discuss the theory relative to experimental data. In order to do this, we need to establish the experimental volumes at zero temperature to enable a proper comparison with the zero-temperature modeling. Thankfully, Söderlind and Young [61] carefully determined these experimental zero-temperature volumes (V 0 ) for Th, Pa, and U (32.55, 24.72, and 20.47 Å 3 , respectively). For neptunium we use the experimental linear coefficient of thermal expansion (27 × 10 −6 K −1 ) [62] combined with the quasi-harmonic Debye-Grüneisen model [4] to back out V 0 (18.9 Å 3 ). For α-plutonium the Debye-Grüneisen procedure is inaccurate because the temperature dependence of the volume is highly anomalous at low temperatures. At temperatures close to 50 K, there is negative thermal expansion observed in dilatation experiments [63]. This anomaly was proposed to be due to magnetism [63], and this speculation is plausible because the negative thermal expansion observed in δ-plutonium is explained from first-principles theory by magnetic interactions [40]. Here, we extract V 0 (19.53 Å 3 ) for α-plutonium from a linear extrapolation, from the three lowest temperatures, to zero temperature of the dilatation data shown in Figure 1. We compare our VASP results, that include spin-orbit interaction and orbital polarization, with the experimental V 0 values in Figure 2. The agreement is very good and almost as good as all-electron (FPLMTO) data that also include spin-orbit coupling and orbital polarization [24]. The small differences are likely due to the PAW method used in VASP.
Appl. Sci. 2019, 9, x FOR PEER REVIEW 6 of 14 of the volume is highly anomalous at low temperatures. At temperatures close to 50 K, there is negative thermal expansion observed in dilatation experiments [63]. This anomaly was proposed to be due to magnetism [63], and this speculation is plausible because the negative thermal expansion observed in δ-plutonium is explained from first-principles theory by magnetic interactions [40]. Here, we extract V0 (19.53 Å 3 ) for α-plutonium from a linear extrapolation, from the three lowest temperatures, to zero temperature of the dilatation data shown in Figure 1. We compare our VASP results, that include spin-orbit interaction and orbital polarization, with the experimental V0 values in Figure 2. The agreement is very good and almost as good as all-electron (FPLMTO) data that also include spin-orbit coupling and orbital polarization [24]. The small differences are likely due to the PAW method used in VASP. Figure 1. Experimental temperature dependence of the atomic volume for α-plutonium. Redrawn from [63] and assuming that the room-temperature (300 K) atomic volume is 20.0 Å 3 . The 0 K value is obtained from a linear extrapolation from the three lowest temperatures (red line).

Figure 1.
Experimental temperature dependence of the atomic volume for α-plutonium. Redrawn from [63] and assuming that the room-temperature (300 K) atomic volume is 20.0 Å 3 . The 0 K value is obtained from a linear extrapolation from the three lowest temperatures (red line). We assess the pseudopotential approximation in Figure 3 where scalar-relativistic VASP and FPLMTO atomic volumes are compiled. Notice that the all-electron FPLMTO volumes are slightly larger, as are the FPLAPW volumes (not shown), but overall the agreement is quite good. For αplutonium the difference is more significant but, in this case, when spin-orbit and orbital polarization are included, the agreement between the methods improves. In our judgement, the VASP pseudopotential is very reasonable for Th to Pu. We assess the pseudopotential approximation in Figure 3 where scalar-relativistic VASP and FPLMTO atomic volumes are compiled. Notice that the all-electron FPLMTO volumes are slightly larger, as are the FPLAPW volumes (not shown), but overall the agreement is quite good. For α-plutonium the difference is more significant but, in this case, when spin-orbit and orbital polarization are included, the agreement between the methods improves. In our judgement, the VASP pseudopotential is very reasonable for Th to Pu.
We thus have a model (DFT + SO + OP), that reproduces the equilibrium volumes for the actinides very accurately. The parabolic behavior of the volumes is duplicated as well as the upturn observed between Np and Pu. The upturn for plutonium is due to several effects: Crystal structure, spin-orbit coupling, orbital polarization, and formation of magnetic moments. The present model does not depend on an ad hoc Hubbard U parameter. Actually, introducing such a parameter in a DFT + U approach invokes physics that is inappropriate for these metals, unless the parameter is very small and inconsequential. To illustrate this, we performed VASP DFT + U calculations that include spin-orbit interaction with the proper exclusion of the p states, as discussed. In Figure 4 we show these results that suggest an over-estimation of all the volumes, except for thorium. Accounting for the fact that the pseudopotential approximation tends to slightly under-estimate the volume (Figure 3), the conclusion from Figure 4 is that the DFT + U volumes are decidedly too large and worse than those of the treatment without a U (Figure 2). The overall effect of the applied U is here relatively small, cf. Figures 2 and 4, because we investigate a very small U (0.4 eV). Using the more common and larger value of U (~1 eV) exacerbates this over-estimation problem for the volumes, while adopting a larger U (>4 eV), often applied for δ-plutonium, causes wholly unrealistic volumes for all the actinides (including δ-plutonium). We thus have a model (DFT + SO + OP), that reproduces the equilibrium volumes for the actinides very accurately. The parabolic behavior of the volumes is duplicated as well as the upturn observed between Np and Pu. The upturn for plutonium is due to several effects: Crystal structure, spin-orbit coupling, orbital polarization, and formation of magnetic moments. The present model does not depend on an ad hoc Hubbard U parameter. Actually, introducing such a parameter in a DFT + U approach invokes physics that is inappropriate for these metals, unless the parameter is very small and inconsequential. To illustrate this, we performed VASP DFT + U calculations that include spin-orbit interaction with the proper exclusion of the p states, as discussed. In Figure 4 we show these results that suggest an over-estimation of all the volumes, except for thorium. Accounting for the fact that the pseudopotential approximation tends to slightly under-estimate the volume ( Figure  3), the conclusion from Figure 4 is that the DFT + U volumes are decidedly too large and worse than those of the treatment without a U (Figure 2). The overall effect of the applied U is here relatively small, cf. Figures 2 and 4, because we investigate a very small U (0.4 eV). Using the more common and larger value of U (~ 1 eV) exacerbates this over-estimation problem for the volumes, while adopting a larger U (>4 eV), often applied for δ-plutonium, causes wholly unrealistic volumes for all the actinides (including δ-plutonium) .   Another very obvious sign of the improper physics that the Hubbard U parameter enforces is found in the prediction of the ground-state phases. For plutonium, DFT + U predicts the wrong ground state, even for the very small U we have adopted in this study (0.4 eV). In Figure 5 we show our fully (both structural and magnetic-configurational) relaxed DFT + U calculations for α-and δplutonium. Notice that the δ-phase energy is about 7 mRy (per atom) lower than the α phase. The fact that the δ phase is predicted to be the ground state when a large U (~ 4 eV) is employed has been discussed [40], but it is remarkable that even a very small U promotes the wrong ground state. Actually, the two phases are energetically degenerate for a small U ~ 0.2 eV. Here, both the α and δ volumes are significantly under-estimated (18.74 and 24.25 Å 3 , respectively). We find that increasing U aggravates the problem and further increases the stability of δ over α (not shown). We emphasize that a theory assuming no ad hoc Coulomb repulsion (U = 0) but orbital polarization and careful Another very obvious sign of the improper physics that the Hubbard U parameter enforces is found in the prediction of the ground-state phases. For plutonium, DFT + U predicts the wrong ground state, even for the very small U we have adopted in this study (0.4 eV). In Figure 5 we show our fully (both structural and magnetic-configurational) relaxed DFT + U calculations for αand δ-plutonium. Notice that the δ-phase energy is about 7 mRy (per atom) lower than the α phase. The fact that the δ phase is predicted to be the ground state when a large U (~4 eV) is employed has been discussed [40], but it is remarkable that even a very small U promotes the wrong ground state. Actually, the two phases are energetically degenerate for a small U~0.2 eV. Here, both the α and δ volumes are significantly under-estimated (18.74 and 24.25 Å 3 , respectively). We find that increasing U aggravates the problem and further increases the stability of δ over α (not shown). We emphasize that a theory assuming no ad hoc Coulomb repulsion (U = 0) but orbital polarization and careful relaxation of all degrees of freedom reproduces the correct ground state in accordance with previous conclusions [18,40].

Summary and Conclusions
We conducted DFT and DFT + U calculations for the purpose of evaluating the potential need for an explicit Hubbard U intra-atomic Coulomb repulsion for the actinide metals Th to Pu. We limit our study to the α phases of these metals (except for δ-plutonium) and find that one carefully needs to address several key approximations in the theory in order to reach a definite conclusion.
First, the spin-orbit interaction has to be introduced in an appropriate fashion. Based on comparisons with Dirac-FPLAPW calculations that include spin-orbit interaction exactly, we recommend a way to accurately account for SO interaction without having to solve the spin-polarized Dirac equation. Namely, by excluding this interaction on the p states, one can accurately apply the usual scalar-relativistic Schrödinger-type equation with perturbative spin-orbit coupling [64] as implemented in VASP and other codes.
Second, one has to make sure that the pseudopotential employed in plane-wave methods is a good representation of the core electrons. We have confirmed that this is the case here by comparisons with results from all-electron codes. The most accurate result is obviously obtained from all-electron approaches, however.
Third, the crystal structures need to be relaxed (i.e., optimized so that the lowest energy is obtained). This is particularly important for α-plutonium because of the sensitive correlation between its crystal and electronic structure. The plutonium computations we have performed are thus nontrivial because forces are required for structural relaxations and these are troublesome to evaluate in

Summary and Conclusions
We conducted DFT and DFT + U calculations for the purpose of evaluating the potential need for an explicit Hubbard U intra-atomic Coulomb repulsion for the actinide metals Th to Pu. We limit our study to the α phases of these metals (except for δ-plutonium) and find that one carefully needs to address several key approximations in the theory in order to reach a definite conclusion.
First, the spin-orbit interaction has to be introduced in an appropriate fashion. Based on comparisons with Dirac-FPLAPW calculations that include spin-orbit interaction exactly, we recommend a way to accurately account for SO interaction without having to solve the spin-polarized Dirac equation. Namely, by excluding this interaction on the p states, one can accurately apply the usual scalar-relativistic Schrödinger-type equation with perturbative spin-orbit coupling [64] as implemented in VASP and other codes.
Second, one has to make sure that the pseudopotential employed in plane-wave methods is a good representation of the core electrons. We have confirmed that this is the case here by comparisons with results from all-electron codes. The most accurate result is obviously obtained from all-electron approaches, however.
Third, the crystal structures need to be relaxed (i.e., optimized so that the lowest energy is obtained). This is particularly important for α-plutonium because of the sensitive correlation between its crystal and electronic structure. The plutonium computations we have performed are thus non-trivial because forces are required for structural relaxations and these are troublesome to evaluate in the presence of strong spin-orbit and orbital-orbital coupling.
The conclusion is that DFT, with the GGA for the electron exchange and correlations, is sufficient for all studied actinide metals (no Hubbard U needed). Recently, numerous others came to the same conclusion [65][66][67][68] for uranium. For plutonium, however, treating the orbital-orbital interaction is necessary for a good electronic structure, total energy, and magnetic form factor [69]. Importantly, we show that DFT + U (U > 0.2 eV) for plutonium predicts the wrong ground-state phase. Because DFT + U defines the static limit of dynamical-mean field theory, that incorporates similar electron correlations, it seems unavoidable that it is also inappropriate for plutonium.
Author Contributions: B.S. and P.S. conceived the study. B.S., P.S., and A.K. performed electronic-structure calculations. All authors helped writing the paper.
Nevertheless, the PAW equations are derived under the assumption of Equation (A4). As a result, the expectation value of any local operatorB confined to the muffin-tin sphere around a nucleus can be simply expressed as follows Ψ kn |B|Ψ kn = ilm,i l m Ψ kn | p ilm φ ilm |B|φ i l m p i l m | Ψ kn . (A5) Finally, when dealing with non-collinear magnetism, the crystal wave functions acquire an extra dimension and become spinor-wave functions Ψ τ kn , with τ now indexing the two spinor components. We are now in a position to derive the expression for the expectation value of any component (α = x, y, z) of the local angular-momentum operator L f α that acts on the f electrons within an atomic sphere. It can be written as where σ ττ α is the Pauli spin matrix in the α direction. The corresponding self-consistent equations are constructed from functional differentiations of Equation (2) with respect to the pseudo spinor wave functions Ψ τ kn and occupations f k,n .