Parameterisation of an LSDA+ model for noncollinear magnetic configurations: Multipolar magnetism in UO
Abstract
Unraveling the complexity of noncollinear magnetism in materials with strongly correlated electrons is a considerable task that requires developing and applying state of the art theories and computational methods. Using the Coury model Hamiltonian, which includes spin and orbital degrees of freedom and generalizes the collinear Stoner Hamiltonian, we derive an extension of the collinear LSDA+ approximation to noncollinear magnetic configurations and explore the magnetic ground state of the archetypal spinorbit correlated oxide UO. We show that parameterizing a noncollinear LSDA+ model requires only one parameter , as opposed to the difference between the Hubbard and Stoner parameters found in an earlier derivation based on a collinear model Hamiltonian. To find the magnetic ground state of UO in the noncollinear configuration space, we combine LSDA+ with a spin adiabatic occupation matrix approach, involving the construction of a magnetic energy surface that follows the adiabatic evolution of the occupation matrix as a function of the spin canting angle. Our results show that the strong spinorbit coupling (SOC) is the key factor stabilizing the socalled 3 spin ordered magnetic ground state of UO. Using a relativistic atomic Hamiltonian we find that the SOC strength is colossal, 1.45 eV per uranium atom, the largest value ever found in relativistic materials. This unusually strong SOC implies that the spin and orbital degrees of freedom are virtually inseparable. As a result, to derive and quantify spinspin interactions it is necessary to adopt the pseudospin picture. By constructing an extended effective multipolar pseudospin Hamiltonian, we prove the significance of octupolar and dipoledipole exchange couplings in establishing the 3 magnetic phase, consistent with the noncollinear spin arrangement, and associated with the canted orbital ordering of uranium orbitals. Importantly, our study reveals that despite the presence of strong spinlattice interaction in UO, the cooperative JahnTeller instability does not contribute to the onset of the noncollinear 3 state, which remains the most favorable ground state even in the undistorted lattice. Finally, we discuss the part played by parameter in the LSDA+ scheme and provide evidence that the choice of this parameter has a substantial quantitative effect on the predicted properties of the oxide, in particular the magnetic exchange interactions and, perhaps trivially, on the band gap: the value of computed fully ab initio by the constrained random phase approximation, =3.46 eV, delivers a band gap of 2.11 eV in good agreement with experiment, and a balanced account of the other relevant energy scales.
I Introduction
A realistic treatment of magnetic properties by ab initio methods requires using approaches that do not artificially constrain the orientation of magnetic moments to a specific direction. This is achieved by means of the socalled noncollinear density functional theory approximation Kubler et al. (1988); Sandratskii (1998). In the noncollinear approximation, the direction of local magnetic moments varies from point to point in real space, and the fact that magnetic noncollinearity does indeed occur in real materials is confirmed by direct experimental observations and first principles calculations Yamada et al. (1970); Heinze et al. (2000); Ohldag et al. (2001); Heide et al. (2009); Hobbs et al. (2000); Hauptmann et al. (2017). Noncollinear spin orderings is particularly evident in heavy element compounds where it results from strong spinorbit coupling effects. Examples include electron systems Faber and Lander (1976); Burlet et al. (1986); Giannozzi and Erd??s (1987) and 5 transition metal oxides Jackeli and Khaliullin (2009); Liu et al. (2015). Fluctuations of magnetic noncollinearity is a factor contributing to electrical and thermal resistivity of alloys Starikov et al. (2018) and may offer a clue to understanding electronic and magnetic phase transitions Kim et al. (2016, 2017a). It has also been discovered that noncollinear magnetic fluctuations may be responsible for the observed anomalous thermal conductivity of uranium dioxide Gofryk et al. (2014), a material of considerable technological significance. Indeed, while in metallic alloys, magnetic fluctuations represent just one of the types of electronic excitations associated with the electric and thermal resistivity of the material, in a semiconducting oxide like UO, where thermal conductivity is dominated by the transport of phonons, the strong spinorbit coupling (SOC) linking atomic displacements with magnetic degrees of freedom provides an additional and potentially significant channel of dissipation. The fact that lattice and magnetic degrees of freedom in uranium dioxide are strongly coupled has been confirmed by observations of piezomagnetism Jaime et al. (2017). The observations were modeled phenomenologically, assuming the occurrence of direct coupling between atomic displacements and magnetic moments. On the electronic structure level, such a coupling stems from the relativistic SOC between the magnetic moments of uranium ions and their orbital degrees of freedom, where the relatively weak but still notable directional character of bonds involving valence orbitals of uranium ions generates interatomic forces that depend on the orientation of magnetic moments.
The fundamental difference between a magnetic metal and an actinide oxide like UO is that while the spinorientationdependent forces in a metal stem from the coordinate dependence of the Heisenberg exchange Ma et al. (2016), in an antiferromagnetic semiconducting actinide oxide the forces result from a combination of Anderson’s superexchange Anderson (1950, 1959), the strong correlations between electrons occupying shells of uranium atoms Kotani and Yamazaki (1992), and the relativistic spinorbit interaction Mironov et al. (2003); Pi et al. (2014a, b) coupling magnetic configurations to forces acting on the nuclei. In what follows, we shall focus on exploring the electronic and magnetic aspects of the problem using a suitably adapted ab initio LSDA+ model.
The subtlety of the problem is associated with the fact that an ab initio treatment of directional magnetic degrees of freedom in materials with strongly correlated electrons requires a noncollinear formulation of the electronic structure model, which at the same time must be suitable for the evaluation of the total energy of the electronic system. A LSDA+ model, often applied to total energy calculations Dudarev et al. (1998), was derived from a collinear model Hamiltonian, similar to the Hamiltonians used earlier by Anisimov et al. Anisimov et al. (1991) and by Kotani and Yamazaki Kotani and Yamazaki (1992). It can be readily shown that these model Hamiltonians are identical to the Hamiltonian of the collinear Stoner model NguyenManh and Dudarev (2009), and hence do not provide a suitable foundation for the treatment of noncollinear magnetism. A LSDA+ model that can in principle be applied to the treatment of noncollinear magnetism was proposed by Liechtenstein et al. Liechtenstein et al. (1995). However the choice of the double counting term, which has the same form as in the collinear case, cf. equation (4) of Ref. Liechtenstein et al. (1995) and equation (3) of Ref. NguyenManh and Dudarev (2009), introduces an element of uncertainty in the total energy part of the analysis.
The choice of the LSDA+ model for the analysis below is stimulated by the extensive application of ab initio approaches to the simulation of materials with strongly correlated electrons Anisimov et al. (1997); Cococcioni (2012); Mazurenko et al. (2007); Ylvisaker et al. (2009); Bousquet and Spaldin (2010); Mellan et al. (2015). In what follows, we explore a generalization of the collinear LSDA+ model Dudarev et al. (1998) to noncollinear magnetic configuration.
The new development has been stimulated and enabled by a recent work by Coury et al. Coury et al. (2016), who found a way of transforming, through an exact calculation, the general fourindex matrix form of the secondquantized Hamiltonian for interacting electrons, into a simpler form, resembling the Hamiltonian used earlier for the derivation of simplified versions of the LSDA+ model Anisimov et al. (1991); Dudarev et al. (1998). The relative simplicity of the Coury Hamiltonian enables performing a more accurate derivation of the LSDA+ model, suitable for the treatment of noncollinear magnetism in materials with strongly correlated electrons, at the same time providing the means for deriving a double counting term for evaluating the total energy.
Below we apply the method to the investigation of electronic structure and noncollinear magnetism in uranium dioxide. This oxide, despite the fact that it has been studied for several decades Allen (1968a, b); Burlet et al. (1986); R. et al. (1999); Laskowski et al. (2004); Zhou and Ozoliņš (2011); Korotin et al. (2016); Pegg et al. (2017); Bruneval et al. (2018), remains a subject of extensive research as the scope of ab initio methods expands to include the treatment of defects Crocombette et al. (2001) and their diffusion Dorado et al. (2012); Vathonne et al. (2017); Liu and Andersson (2018). It is known that in metals the structure of defects NguyenManh et al. (2006) and atomic diffusion Wen et al. (2013) are highly sensitive to magnetism. The experimental observation showing the effect of magnetic fluctuations on the thermal conductivity of UO Gofryk et al. (2014) suggests that other hightemperature transport properties may also be affected by the magnetic degrees of freedom.
The magnetic ground state of defectfree UO has been extensively investigated both experimentally and theoretically Allen (1968a, b); Burlet et al. (1986); R. et al. (1999); Laskowski et al. (2004); Zhou and Ozoliņš (2011); Caciuffo et al. (2011). However, a definitive verdict on the stability of the 3 magnetic ground state structure and the form of a lowenergy spin excitation Hamiltonian suitable, for example, for spinlattice dynamics simulations Ma et al. (2008, 2016); Tranchida et al. (2018), still remains outstanding. Partly, the difficulty is associated with the fact that an ab initio treatment of the oxide requires exploring the complexity of the energy landscape involving multiple local energy minima, which impedes applications of conventional energy minimization algorithms Dorado et al. (2013). Also, it is necessary to take into account relativistic effects, giving rise to large SOC Korotin et al. (2016) and the emergence of multipolar spin interactions. Furthermore, it is necessary to determine the population of uranium orbitals, which result in the formation of orbital magnetic moments on uranium ions that are twice the magnitude of the spin moment Dudarev et al. (2000). Finally, magnetism is intimately linked to distortions of the underlying crystal structure, manifested in UO via the JahnTeller (JT) instabilities Sasaki and Obata (1970), and noncollinear spin arrangement in the form of a 3ordered magnetic state. In this paper, we address the outstanding issues noted above using a consistent theoretical and methodological framework, aiming to uncover the ultimate quantum origin of the magnetic ground state in UO.
An important aspect of application of LSDA+ models to ab initio simulations of materials is the proper choice of values of the Coulomb repulsion and Hund’s coupling parameters and . Owing to the difficulties inherent in evaluating the magnitude of these interactions, these quantities are generally treated as tunable parameters, and are chosen according to how well they reproduce some specific properties such as the band gap, the equilibrium volume, the magnetic moment or the formation energy Franchini et al. (2007). In an attempt to improve the predictive capacity of LSDA+ models, several approaches have been proposed that aim to compute and without relying on phenomenological considerations Aryasetiawan et al. (2006); Cococcioni and de Gironcoli (2005); Solovyev (2006); Anisimov and Gunnarsson (1991). Among them, the constrained random phase approximation (cRPA) allows one to take into account the effect of screening of the Coulomb interaction between correlated electrons and to estimate the magnitude of parameters and in the correlated subspace. The cRPA method has been successfully applied to a large variety of systems, and now represents a reliable and computationally viable method of obtaining accurate values of and Vaugier et al. (2012); Kim et al. (2017b, 2018); Liu et al. (2018). In the analysis below, we apply this approach to computing and in the correlated subspace, in this way providing a solid reference for the choice of values of these parameters in actinide compounds, and allowing for a fully ab initiobased application of the LSDA+ model. In particular, the possibility to selectively compute and enables making an estimate of the role played by in the parametrization of the LSDA+ model for noncollinear spin systems.
The manuscript is organized as follows. We start from deriving a noncollinear magnetic LSDA+ model, giving equations for the effective correcting oneelectron potential and the double counting term, showing that the relevant formulae require only one parameter as opposed to the difference found in the collinear version of the model. We then use cRPA to compute the value of and compare the results with the available data. Subsequently, we address the question of stability of the noncollinear 3 magnetic structure and, by using an adiabatic occupation matrix approach, prove that the 3 structure represents the lowest energy configuration of UO even in a perfect cubic lattice. Finally, we perform a detailed analysis of effective spin Hamiltonians and spinspin interactions, concluding with the assessment of potential applications of the proposed methodology to computing high temperature properties, including simulations of atomic and magnetic dynamics.
Ii a Noncollinear LSDA+ model
The existing collinear LSDA+ model, which aims to provide an improved description of the electronic structure of materials characterized by strong electron correlations in spatially localized d and f shells, adds a correction term to the effective single particle electron potential Anisimov et al. (1991); Liechtenstein et al. (1995); Dudarev et al. (1998)
(1) 
and a double counting correction to the total electronic energy Dudarev et al. (1998)
(2) 
The latter is required to take into account the fact that a simple sum of single particle energies in a system of interacting electrons does not represent the correct contribution to the total energy, see e.g., Eq. (15) and Eq. (16) of Ref. Dudarev and Derlet (2007). In Eq. (1) and Eq. (2), indexes refer to orbitals associated with a lattice site, and is the spin index.
Eq. (1) and Eq. (2) were derived from a model tightbinding Hamiltonian Caroli et al. (1969); Kotani and Yamazaki (1992), where the onsite electron interaction terms have the form
(3) 
It can be shown NguyenManh and Dudarev (2009); Coury et al. (2016) that the above Hamiltonian is identical to the Hamiltonian of the collinear Stoner model
(4) 
where , and . Using the procedure outlined in Ref. Dudarev et al. (1998), Eq. (1) and Eq. (2) can be derived by evaluating the expectation values of either Eq. (3) or Eq. (4).
Despite the relative success of the LSDA+ model Cococcioni (2012), there are two points in the derivation that remain elusive. It is unclear to what extent the approximations associated with Hamiltonian (3) affect the form of Eq. (1) and Eq. (2), and also how to generalize these equations to noncollinear magnetic configurations, significant to applications Zhou and Ozoliņš (2011); Gofryk et al. (2014); Ma et al. (2016). The brief analysis given below addresses these two points.
First, we note that there is a major term missing in collinear Hamiltonians (3) and (4). This missing term has been identified in Ref. Coury et al. (2016), and we show below that this term contributes to the LSDA+ correction, changing its form. Second, we note that the correction itself is not invariant in the extended space of spin and orbital indexes, a point that can be readily rectified using a suitable definition of the convolution of the full spin and orbitaldependent electron density matrix.
To derive a noncollinear form of the LSDA+ model, we start by considering a model Hamiltonian, different from and more general than Eq. (3) or Eq. (4), and apply it to the parametrization of the LSDA+ equations, simplifying the full matrix representation of electron interactions. The treatment is based on a recent study by Coury et al. detailed in Ref. Coury et al. (2016). We find that application of the same methodology as what was followed in Ref. Dudarev et al. (1998), to a more complete Hamiltonian Coury et al. (2016), results in the LSDA+ equations that require only one parameter as opposed to the two parameters and entering the formula derived in Dudarev et al. (1998). The general functional form of the LSDA+ correction remains the same.
An onsite Hamiltonian, describing interaction between electrons occupying orbitals , is given by a sum of combinations of four creation and annihilation operators multiplied by a fourindex matrix :
(5) 
Matrix has (2l+1) elements, which in the case of p (l=1) electrons amounts to =81 elements, and =625 elements in the case of d electrons. Symmetry constraints show that all the elements of can be parameterized using only two independent constants in the pelectron case, three constants in the delectron case, and four in the felectron case. In the pelectron cubic harmonic orbital case, using an analogy with the theory of isotropic elasticity Kosevich (1999), where the fourindex matrix of elastic constants has the same symmetry as , Hamiltonian (5) can be written exactly as Coury et al. (2016)
(6) 
Here, is the operator of the total number of electrons on a site, , and
is the total vector magnetic moment operator associated with the site. In Eq. (6), :: denotes normal ordering of creation and annihilation operators. The normally ordered terms in Eq. (6), expressed using conventional notations, have the form
(7) 
The first two terms in Eq. (6) are identical to those entering Eq. (4), with the natural exception that now the magnetic moment operator is a vector quantity. Most significantly, in addition to the terms contained in Eq. (4), Hamiltonian (6) includes an extra third term, required by symmetry and absent in Eq. (4). This term is related to the orbital moment of electrons on a lattice site Coury et al. (2016). In the collinear approximation, where operator in Eq. (6) is replaced by its projection on a quantization axis, and the third term in Eq. (6) is neglected, Hamiltonian (6) is identical to Eq. (4). We note that Hamiltonian (6) is exact in the sense that no procedure of “directional” averaging is involved in the transformation from Eq. (5) to Eq. (6).
We now follow the derivation given in Refs. Anisimov et al. (1991); Dudarev et al. (1998) and deduce the LSDA+ model from Hamiltonian (6). We identify the terms in Hamiltonian (6) that contain two creation and two annihilation operators acting on the same electronic state . In the meanfield approximation, these terms provide contribution to the total energy proportional to whereas their exact expectation value is proportional to . The LSDA+ model correction equals the difference between the exact and meanfield expectation values of these terms, resulting in
(8) 
In the above expression, each term in square brackets corresponds to a respective term in Hamiltonian (6), and is the electron occupation number of an orbital state with spin index .
We see that the terms in Eq. (8), which contain parameter , cancel each other, and only the term proportional to parameter remains. The sum of the first two terms in square brackets gives the coefficient found in the LSDA+ derivation given earlier Dudarev et al. (1998), which was based on a collinear model Hamiltonian (4). The third term in square brackets in Eq. (8) stems from the last term in Hamiltonian (6), which is missing in the Hamiltonians Caroli et al. (1969); Kotani and Yamazaki (1992) used earlier. The full cancellation of the terms containing parameter is perhaps not surprising, given the original form of the Hubbard Hamiltonian Hubbard (1963) that contains no terms and is still able to generate a variety of magnetic solutions, originating solely from strong onsite electron correlations.
A general form of Eq. (8), invariant with respect to the choice of electronic orbitals and spin quantization axis, is now:
(9)  
where is the full orbital and spindependent oneelectron density matrix.
In the collinear approximation, where the density matrix is assumed to be diagonal with respect to the subset of its spin indexes , Eq. (9) is identical to equation (5) of Ref. Dudarev et al. (1998). In this limit, the only difference between equation (5) of Ref. Dudarev et al. (1998) and Eq. (9) given above is the prefactor that, due to the presence of the third term in Hamiltonian (6), now equals as opposed to derived from Hamiltonian (3).
An invariant orbital and spindependent noncollinear form of LSDA+ (1) and (2) is now:
(10) 
and
(11) 
The need to add the double counting term (11) to the total energy, evaluated using the conventional KohnSham procedure, stems from the fact that the singleparticle electron potential (10) depends on the occupancy of electron orbitals through the term proportional to . It is this occupancydependent term in Eq. (10) that makes a simple sum of oneparticle energies different from the total energy given by Eq. (9). The above equations show that, in addition to correcting the prefactor in the formula, a fully invariant LSDA+ construction requires convoluting the density matrix over the full set of its orbital and spin indexes, a point that was not addressed in earlier derivations Anisimov et al. (1991); Liechtenstein et al. (1995); Dudarev et al. (1998).
Similarly to Eq. (8), the terms containing parameter also cancel exactly if we perform the above derivation for the delectron case Coury et al. (2016). The most direct way of carrying out the derivation is to start from equation (22) of Ref. Coury et al. (2016) and note that in the delectron case, all the terms containing parameter can be expressed in terms of a renormalised parameter , resulting in Eq. (9) given above, plus small terms proportional to , which together amount to a small fraction of an electronvolt per atom and are normally neglected in applications Cococcioni (2012). This suggests that the singleparameter form of the LSDA+ correction given by Eq. (10) and Eq. (11) remains sufficiently accurate and applicable to delectron orbitals, and other types of shells containing correlated electrons. Parameter , according to the analysis given in Anisimov et al. (1997); Vaugier et al. (2012), is an effective quantity, characterizing the strength of electronelectron interactions, modified by effects of manybody selfscreening. The significance of the above equations, in addition to the fact that they have been derived from a more accurate form of the onsite interaction Hamiltonian, is in that they can be used for treating noncollinear magnetic configurations.
Iii Ab Initio Methodology
This section details the set of theoretical and methodological approaches adopted in this study. All the calculations were carried out using the Vienna ab initio simulation package (VASP) Kresse and Hafner (1993); Kresse and Furthmüller (1996), where the noncollinear extension of the Dudarev et al. LSDA+ scheme is implemented Bengone et al. (2000) with the full inclusion of relativistic effects (e.g., a full selfconsistent treatment of spinorbit coupling) Hobbs et al. (2000). A robust energy cut off up to 700 eV with the convergence precision of 10 eV was used in all the calculations, and the Brillouin zone was sampled using a 666 point mesh. Atomic positions were optimized with the lattice parameters fixed at its observed value (=5.469 Å) Momin et al. (1991). In what follows, we explore in considerable detail (A) the evaluation of interaction parameters and within the cRPA, (B) magnetically constrained density functional theory (DFT), (C) spin adiabatic occupation matrix approach to the calculation of magnetic energy landscape, (D) effective multipolar pseudospin Hamiltonian and exchange coupling, and (E) a strategy for estimating the strength of the SOC.
iii.1 Constrained random phase approximation
Interaction parameters and were computed from first principles using the constrained random phase approximation Aryasetiawan et al. (2006). In the cRPA, the Coulomb repulsion and Hund’s coupling parameters and are derived from the matrix elements of written in terms of Wannier basis functions representing the correlated subspace (in this case, the uranium states)
(12) 
More precisely, and are the average matrix elements (intra and interorbital interactions of the different spins) and (interorbital interactions of the same spins), respectively. In Eq. (12), is the partially screened interaction kernel, which is calculated by solving the Dysonlike equation
(13) 
where is the bare (unscreened) interaction kernel and is the polarizability, excluding contributions from the “target” correlated subspace, . Following this procedure, we have found =3.46 eV and =0.30 eV, corresponding to the effective interaction parameter =3.16 eV. These values are expectedly smaller than those extracted from experimental estimates, adopted in previous LSDA+ studies of UO Dudarev et al. (1997): =4.50 eV and =0.54 eV, =3.96 eV. In order to verify the role played by the prefactor entering Eq. (9) in the description of magnetic properties of UO, we have tested four different choices of interaction parameters, namely

=3.16 eV (i.e. )

=3.46 eV (i.e. )

=3.96 eV (i.e. ‘Expt.’ )

=4.50 eV (i.e. ‘Expt.’ )
iii.2 Magnetically constrained noncollinear DFT+
To determine the noncollinear magnetic ground state of UO, we have minimized the total energy, treating it as a function of the direction of spin moments and using the magnetically constrained noncollinear DFT+ Hobbs et al. (2000); Liu et al. (2015); Ma and Dudarev (2015). In particular, we have inspected spin rotations that transform the system from a characteristic noncollinear 3 state into collinear antiferromagnetically (AFM) ordered and configurations Zhou and Ozoliņš (2011), illustrated in Fig. 1. A noncollinear 3 phase is described by three independent wave vectors and can be represented by a combination of three different phases, one longitudinal and two equivalent transverse. To facilitate the construction of the canted magnetic energy landscape, we used the longitudinal 3 ordered magnetic structure shown in Fig. 1 as a starting configuration. The two other ordered AFM configurations, and , belong to the 1k (one wavevector) and 2k (two wavevectors) categories, respectively.
The –3– magnetic structure transformation pathway can be defined by a concerted variation of angle on the four inequivalent uranium sites in a UO magnetic unit cell [see Fig. 2(a) and (c)]. Constrained energy minimization as a function of along the transformation pathway is achieved by considering the energy penalty arising from constraining the direction of the spin magnetic moment, defined by the formula
(14) 
Here is the unconstrained DFT total energy, whereas the second term is a penalty contribution defined as a noncollinear directional constraint on the direction of local moments with respect to an arbitrary set of unit vectors on sites . is the magnetic moment computed by integrating over a WignerSeitz cell centered on atom (the effective WignerSeitz radius is 1.588 Å for a U ion and 0.82 Å for an O ion). Parameter defines the magnitude of the energy penalty term, where by progressively increasing , functional (14) is driven to convergence towards the DFT total energy Ma and Dudarev (2015). We used the value =10 that guarantees that the expectation value of the energy penalty term is lower than 10 eV.
iii.3 Adiabatic spin occupation matrix approach
A known drawback of DFT+ approaches is the difficulty associated with finding the lowest energy state of a strongly correlated magnetic material. In most cases a DFT+ functional exhibits a multitude of local minima corresponding to a variety of spin and orbital occupancies in the correlated electronic subspace Dorado et al. (2009); Meredig et al. (2010); Zhou and Ozoliņš (2011); Allen and Watson (2014). The difficulty with finding a global minimum stems from the curvature of the energy surface as a function of orbital occupations Meredig et al. (2010). In standard DFT calculations the energy surface is typically convex, but the global minimum might correspond to a physically unreasonable partial fractional orbital occupation, e.g. a metallic state of the material that in reality is an insulator. DFT+ corrects this by adding a term that penalizes fractional occupations, but this correction changes the curvature of the energy surface from convex to concave, unavoidably giving rise to many local minima Meredig et al. (2010). Dorado and coworkers proposed to address this point by performing a search involving a large number of selfconsistent calculations, each starting from different initial occupation matrices, and to select the outcome corresponding to the lowest total energy Dorado et al. (2009). This procedure can be accelerated by adiabatically “turning on” the value of parameter starting from the DFT limit and gradually converging to the true ground state with integral orbital occupations Meredig et al. (2010). The above issue is particularly pertinent to noncollinear spin systems, where small rotations of spin moments could give rise to many local minima, all lying in a few meV energy range.
Bearing in mind this aspect of energy minimization, we have combined the occupation matrix approach with a gradual adiabatic change of the spin moment direction, using the magnetically constrained noncollinear DFT+ functional described in the previous subsection Hobbs et al. (2000); Liu et al. (2015); Ma and Dudarev (2015). Starting from the 3type noncollinear spin ordered state of UO Burlet et al. (1986); Zhou and Ozoliņš (2011) shown in Fig. 2(a), we gradually change the canting angle , moving adiabatically from the noncollinear 3 state to the energetically comparable AFM collinearly ordered and configurations shown in Fig. 2(c). At each canting step, we initialize the occupation matrix to the one obtained at the preceding step and perform a fully selfconsistent calculation. In this way, by gradually perturbing the wave function of the 3 state, we are able to construct a smooth total energy curve E() as a function of the canting angle , shown in Fig. 2(c). The absence of cups and sudden jumps guarantees that this energy curve represents the lowest energy path linking the spin configurations considered here, and also that the 3 state is indeed the global minimum with respect to spin rotations. Probing other spin configurations (not shown here) confirmed this analysis and did result in finding spin configurations with lower total energies. We note that the orbital moment remains antiparallel to the spin moment everywhere on the transformation pathway, and the magnitudes of both moments are almost independent of the canting angle: and [see Fig. 2(b) and (c)].
iii.4 Effective pseudospin Hamiltonian and exchange interactions
To characterize magnetic properties of a material and to understand the origin of the specific spin ordered configuration that it adopts, it is necessary to quantify the dominant spinspin interactions. For systems conserving the total spin moment, magnetic coupling parameters can be analyzed in terms of an effective Heisenberg spin Hamiltonian involving conventional spin operators. In materials with strong spinorbit coupling, the spin moments alone are not conserved and it is more appropriate to use the pseudospin operators and pseudospin Hamiltonian Arima et al. (1969) suitable for the treatment of multipolar interactions Pryce (1950); A. Abragam (1951); Chibotaru (2013),
In the pseudospin picture, spin and orbital degrees of freedom are not independent, and the operator set is formed by a unit multipole (tensor) operators Sakurai (1993); Santini et al. (2009); Pi et al. (2014b) (here is the angular moment, the rank and ). In a general form, a multipolar exchange Hamiltonian can be written as Santini et al. (2009); Pi et al. (2014b)
(15) 
where , are the site indexes, and are the coupling constants describing how the energy of the system changes as a result of variation of the two multipole moments and .
UO adopts a noncollinear 3ordered magnetic configuration, and the twoelectrons () ground state of a uranium ion is a triplet, corresponding to the effective spin (pseudospin) Mironov et al. (2003). The ground state is associated with cooperative quadrupolar interactions that cannot be accounted for by using a Heisenberg model Pi et al. (2014b); Giannozzi and Erd??s (1987); Mironov et al. (2003); Carretta et al. (2010); Caciuffo et al. (2011), but can be modeled by means of a suitable pseudospin Hamiltonian.
To compute the spinspin interactions in UO, we adopt the multipolar spin Hamiltonian derived by Mironov et al. Mironov et al. (2003), describing superexchange interactions between neighboring U ions in the 5 configuration. The general form of the Mironov exchange Hamiltonian reads
(16) 
where is a spinindependent parameter, whereas the remaining terms accounts for various types of spin interactions, which can be written using the conventional spin variables as
(17)  
(18)  
(19)  
(20) 
Here, is a singlespin term, quadratic in the spin components and accounting for the zero field splitting (ZFS) dipolar interactions; the ZFS parameters and describe the axial and transversal components of magnetic dipoledipole (DD) interaction, respectively. is bilinear in spins (and isotropic under rotations) and describes spin exchange interactions (sometimes called the Heisenberg exchange), parameterized by , and . The term describes fourspin exchange interactions with as the corresponding coupling constant. Finally, accounts for biquadratic quadrupolequadrupole (QQ) interactions, where are the components of the quadrupole operator, specifically:
(21)  
(22)  
(23) 
where A, B. Labels A and B refer to the interacting nearest neighbour uranium ions, see Fig. 3 for details. There are four inequivalent sites in a magnetic unit cell of UO [U1–U4, see Fig. 2(a)] producing six inequivalent nearest neighbor AB pairs: A=U1, B=U2; A=U1, B=U3; A=U1, B=U4; A=U2, B=U3; A=U2, B=U4; A=U3, B=U4. and are the two spins forming a distinct inequivalent AB pair, and , , and are the local quantization axes illustrated in Fig. 3. Cartesian components of spins and in the local quantization axis representation are given in the Supplementary Material SM ().
A fully ab initio evaluation of the above superexchange parameters is a difficult task Pi et al. (2014b, a). We also note that depending on the specific definition of tensor operators, slightly different forms of superexchange coupling can be found in literature Mironov et al. (2003); Pi et al. (2014b); Giannozzi and Erd??s (1987); Carretta et al. (2010), impeding an accurate quantitative comparison with existing data.
In formulating the effective pseudospin Hamiltonian we followed the procedure by Mironov Mironov et al. (2003). To estimate the magnitude of effective magnetic interactions, Mironov et al. used a secondorder perturbation approach, using freeion and cubic crystalfield parameters and limiting the interaction to the two nearestneighbor U atoms only. In our analysis, we estimate the dominant exchange couplings by means of a controlled fitting procedure, involving the mapping of onto DFT++SOC total energies using Eq. (16). To achieve this, we have rewritten the four terms entering Eq. (16) as functions of the canting angle , replacing components of spins by their explicit expressions in terms of local Cartesian components and arriving at the total magnetic energy expressed as a function of . After some algebra, we find that has the form
(24) 
where the coefficients are given in terms of the eleven superexchange parameters entering the Mironov Hamiltonian (, , , , , , , , , , and ):
(25)  
(26)  
(27)  
(28)  
(29)  
(30)  
(31)  
(32) 
iii.5 Calculation of the spinorbit coupling strength
We conclude the methodological section with an estimate of the strength of SOC in UO. To produce this estimate, we relate the relativistic total energies obtained from first principles calculations to the relativistic atomic Hamiltonian for orbitals:
(33) 
where defines the SOC strength. By using the 14 () spinors as a basis in the following order: , , , , , , (and the corresponding spinors), the atomic Hamiltonian can be written as a (1414) matrix, see Eq. (34) below Jones and Albers (2009).
The diagonalization of matrix (34) yields the following eigenvalues , , , , , , , , , , , , , . From the eigenvalues of the SOC Hamiltonian, it is possible to extract the contribution to the total energy from the SOC effects, , by considering either the SOCinduced splitting () or the energy contribution arising from the occupied states. We have followed the latter route, as in this case a suitable mapping can be constructed between the atomic limit () and realistic ab initio calculations. Specifically, considering that U ions in UO are in the electronic configuration, the two electrons occupy the lowest two eigenvalues () resulting in . An estimate of can be obtained from the DFT total energy difference between a relativistic (with SOC) and a nonrelativistic calculation (no SOC), i.e., . To exclude the spurious energy contributions arising from differences in the electronic ground states (insulating vs. metallic), this estimate was obtained using and the 3k spin ordered configuration. In this limit =0, both reference states (with and without SOC) are metallic. On the other hand, a DFT++SOC calculation delivers an insulating solution, where DFT+ (without SOC) stabilizies an insulating ground state: this would include terms other than SOC to the energy balance affecting a proper valuation of the SOC energy. Using the above approach, we find that eV per uranium ion, corresponding to =1.45 eV. We should mention that in the metallic solution, the values of spin and orbital moments are greatly reduced with respect to those found in the DFT++SOC ground state (referred to in Table 1), specifically 1 and .
(34) 
Further support for this large value of comes from an approximate scaling of the magnitude of SOC at atomic level, where it is known that the SOC parameter is proportional to , where is the atomic number. By rescaling the SOC strength of iridium (0.5 eV) Andlauer et al. (1976) with the relative nuclear charge of Ir () and U () we find
confirming the very large expected value of the SOC parameter in uranium dioxide.
Iv Results and Discussion
The large SOC in UO is responsible for the formation of a triplet described by an effective pseudospin Mironov et al. (2003) state, where the spin and orbital moments are ordered in a 3 magnetic structure Burlet et al. (1986); Amoretti et al. (1989), see Fig. 2(a) and (b). Moreover, the various types of (multipolar) superexchange interactions acting in the 3 magnetic configuration are coupled with the cooperative JahnTeller (JT) effect, manifested by a distortion of the oxygen cage around the U ions Allen (1968b); Santini et al. (2009); Giannozzi and Erd??s (1987); Jaime et al. (2017). The computational verification of these experimental observations and their interpretation on a quantum level is a difficult task due to a variety of factors: (i) magnetic noncollinearity, (ii) selfinteraction acting in the U manifold, and (iii) existence of multiple local minima in a narrow energy interval Laskowski et al. (2004); Dorado et al. (2009); Zhou and Ozoliņš (2011). As was noted above, a combination of fully relativistic and magnetically constrained DFT+ with the adiabatic evolution of the occupation matrix is able to predict the ground state of UO [Fig. 2(c)] and should help decipher the subtleties of electronic and magnetic effects in UO.
To gain insight into the nature of magnetic interactions, we compute the magnetic energy curves similar to the one shown in Fig. 2, but this time performing the calculations for several different values of parameter used in the DFT+ formalism. The curves shown in Fig. 4 suggest that the noncollinear 3 ordering remains the lowest energy state for any value of , but the energy difference between the 3 phase and the competing AFM collinear phases and , illustrated in Fig. 1, depends sensitively on the choice of . As the value of parameter increases, the relative stability of the 3 state decreases, and it becomes progressively less energetically costly to rotate the spins. This implies that the value of the magnetic exchange interactions is also sensitive to the choice of . We shall discuss this later in the section.
=3.16  =3.46  =3.96  =4.50  

1.52  1.53  1.54  1.54  
3.18  3.19  3.21  3.22  
1.66  1.66  1.67  1.68  
1.91  2.11  2.44  2.78 
The role of parameter is also reflected in the fundamental electronic and magnetic properties of the 3 ground state. Table 1 gives values of the spin moment , the orbital moment , and the insulating gap computed for =3.16, 3.46, 3.96 and 4.50 eV. While the moments are only marginally affected by the choice of (and all of them compare well with the observed total moment of 1.74 Frazer et al. (1965)), the band gap varies significantly from 1.91 eV to 2.78 eV. The best agreement with experiment (2.0 eV Schoenes (1980) ) is found for the relatively small , in agreement with the first principles estimate of the Coulomb interaction parameters based on cRPA (see Sec. III.1), and also in agreement with results derived from a recent fitting analysis Pegg et al. (2017). We also note that even though the value of parameter in UO is not very large, 0.30 eV, reducing by 0.3 eV reduces the band gap by about 10% (=1.91 eV for =3.16 eV and =2.11 eV for =3.46 eV, both fairly close to the observed band gap of approximately 2.0 eV).
Next, we examine the quantum mechanism responsible for the onset of magnetic 3 ordering. We remind the reader that in a JahnTelleractive material with strong SOC, the JahnTeller (JT) instability and exchange interactions are antagonists since the JT effect tends to stabilize states with quenched orbital momentum whereas SOC tends to maximize the orbital momentum Goodenough (1963). This conclusion is generally valid if the crystal field is large, but in UO the strength of SOC is very large ( 1.5 eV according to the estimate above) exceeding the energy scale of crystalfield excitations (150180 meV Amoretti et al. (1989)), and therefore the spinorbit interaction can be safely considered as the dominant energy scale and the leading factor stabilizing the 3 state. To verify this hypothesis, we have calculated the total energy of 1, 2, and 3 magnetically ordered states as a function of strength of the JT effect, switching it on and off. In the on mode we have tested both the selfconsistently derived (0.003 Å) and the experimentally observed (0.014 Å Faber et al. (1975)) values of JTdisplacements, and examined two values of SOC, the full SOC strength eV and half SOC strength eV. The data, given in Table 2, show that the JT effect has virtually no effect on the relative energies of magnetic configurations. The energy landscape and the energy difference between the 3, 2k and 1k states are not affected by the strength of JT distortion and remain essentially unchanged (see Table 2 and Supplementary Material SM ()). This unequivocally demonstrates that the JT effect is not the mechanism that drives the system towards the 3 ground state Pourovskii and Khmelevskyi (2018). On the other hand, rescaling the SOC strength to half the original value ( eV) causes a huge energy change, favorable for both the AFM and AFM configurations, where the latter as a result of halving the SOC value becomes the most favorable one by more than 11 meV/f.u. This provides a clear indication that SOC is the major driving force responsible for the stabilization of the 3 state; reducing the SOC strength leads to the overstabilization of collinear magnetic structures.
JT effect  

3.86  1.01  
3.87  1.03  
3.88  0.94  
SOC strength  
3.87  1.03  
0.5  8.78  11.28 
To discern how SOC stabilizes the 3 ordering of moments, we have explored the differences between electronic and magnetic properties of these three phases. Surprisingly, there are only marginal changes in the magnitude of spin and orbital moments (see Table 3), in the band structure (see Fig. 5), as well as in the occupation numbers of states in the manifold (see Table 4). A closer inspection of the band structure shows that even though the overall bonding picture in all the three magnetically ordered configurations is almost identical (including the size of the band gap), the manifold in the 3 phase exhibits larger SOCinduced splitting, which causes a change in the band topology.
AFM  AFM  3  

1.55  1.51  1.53  
3.26  3.24  3.19  
1.70  1.73  1.66 
To better understand the significance of this change in the topology of band structure, in Fig. 6 we plot chargedensity isosurfaces for the occupied orbitals in the energy interval (2, 0) eV, projected onto a (110) plane containing both uranium and oxygen ions. The results show that the 3 configuration is the only spin arrangement exhibiting a visible orbital anisotropy at the uranium sites, associated with the canted ordering of orbitals. The orbitals are rotated towards the nearest oxygen sites, following the same chessboard configuration of the 3 spin ordering as that shown in Fig. 2(a) and (b). Remarkably, the effect of SOC, critical to the stabilization of the 3 phase, is manifested primarily in the shape of orbitals rather than in the total orbital occupation, as illustrated in Table 4. The energy required to stabilize the 3 state over the AFM and AFM configurations is gained from a SOCinduced rotation of occupancies of particular orbitals, which follows the rotation of the local spin moments and enables constructive interaction with the oxygen electronic states.
AFM  3  AFM 

0.0268  0.0271  0.0273 
0.0274  0.0282  0.0275 
0.0296  0.0288  0.0296 
0.0316  0.0341  0.0333 
0.0348  0.0356  0.0359 
0.0390  0.0365  0.0362 
0.0397  0.0366  0.0375 
0.0398  0.0384  0.0391 
0.0454  0.0488  0.0477 
0.0501  0.0508  0.0513 
0.1233  0.1238  0.1238 
0.1393  0.1407  0.1404 
0.9852  0.9846  0.9846 
0.9888  0.9858  0.9860 
Having established the DFT++SOC as a suitable theory for the ground state electronic properties of UO, highlighting the significance of SOC in this compound, we are now ready to proceed to the analysis of superexchange spin interaction mechanisms, to deduce information about the quantum origin of the 3 state. As we noted in the computational section, this can be done by fitting the magnetic energy computed using ab initio methods, to the multipolar Hamiltonian (24). The presence of eleven parameters in the multipolar Hamiltonian clearly poses a well known problem for the multiparameter fitting procedure von Neumann (). To handle this complication, we rely on the analysis by Mironov et al. Mironov et al. (2003) that can be summarized as follows. First, we note that using the values of parameters evaluated by Mironov (collected in Table 5) in the pseudospin Hamiltonian, already leads to an overall fairly good account of the first principles magnetic energy, as illustrated graphically in Fig. 7. Even though the two curves do not match well, Mironov’s parameters predict the correct position of the energy minimum, located at the 3 position, suggesting that all the relevant magnetic coupling terms are correctly included in the theoretical treatment. However, Mironov’s parameters deliver a curve that varies over a significantly narrower energy interval than the curve derived from the first principles data, and as a result the relative stability of the 3 state with respect to the collinear 1k and 2k states is underestimated by about 50%.
The exchange parameters in Mironov’s model fall into four different categories, namely (i) singlespin parameters and , (ii) bilinear parameters , and , (iii) parameters describing the fourspin terms and and (iv) parameters of biquadratic interactions , , , and . The accuracy of Mironov’s approach can be improved by noting that, according to the calculations by Savrasov and coworkers, the strength of quadrupolar (QQ) interactions computed by Mironov is underestimated by an order of magnitude Pi et al. (2014b). Following this argument, we have fitted by varying the two largest quadrupolar terms (, ) only, and keeping all the other superexchange parameters fixed to the original Mironov’s values. The resulting curve is in excellent agreement with first principles energies (=0.9989, see Fig. 7), and this improvement is associated with a very large increase of the magnitude of the quadrupolar terms, approximately by an order of magnitude (+)/2=4.56 meV, see Table 5. However, this is in very good agreement with the earlier DFT+ data (3.1 meV Pi et al. (2014b)) and the values extracted from experimental spinwave spectra (1.9 meV) Caciuffo et al. (2011). As expected, the values of and are sensitive to the choice of parameter , see Table 6. The strength of the QQ interaction increases as a function of , this in particular applies to the the anisotropic biquadratic interaction . Nevertheless, the resulting values do not depend on the JahnTeller distortions; using the undistorted cubic phase one obtains essentially the the same values of parameters, further demonstrating the fairly negligible role played by the JT effect in stabilizing the 3 noncollinear state.
Mironov Mironov et al. (2003)  0.57  0.64  1.82  2.74  2.33  0.04  0.05  0.22  0.32  0.07  0.02  

Fit (QQ)  7.18  1.94  0.9989 
(+)/2  

= 3.16 eV  6.55  0.36  3.46 
= 3.46 eV  7.18  1.94  4.56 
= 3.96 eV  7.64  3.77  5.71 
= 4.50 eV  7.81  5.19  6.50 
We conclude the discussion of magnetic properties of UO by analyzing the individual contributions of various types of superexchange mechanisms to the stabilization of the 3 state. As was noted in the section on computational methods, the total magnetic Hamiltonian is expressed as a sum of four terms, each corresponding to a specific type of superexchange interaction: accounts for the DD interaction [Eq. (17)], describes the bilinear exchange [Eq. (18)], represents the fourspin exchange [Eq. (19)], and finally takes care of the quadrupolar coupling [Eq. (20)]. The data, summarized in Fig. 8, clearly show that the formation of the 3 state occurs as a result of a concerted action of the DD and octupolar interactions. The contributions of bilinear exchanges (, ’s) and four spin exchanges (, ’s) are essentially independent of the canting angle, resulting in the rather flat curves. On the other hand, the DD interactions have a quadraticlike trend with a marked minimum at 3 and the fitcorrected quadrupolar term (the right panel of Fig. 8) is not only minimum at 3, but also correctly describes the energy pathway from the 3 to the 1 and 2 states, following the trend exhibited by the first principles energies (see Fig. 7).
Based on the above results, we can conclude that the onset of the 3 state in UO is driven, at the quantum level, by a concerted action of the DD and QQ spin interactions. These interactions are active in the undistorted and JTdistorted crystal lattices, clearly indicating that, despite the existing coupling between the spin and lattice degrees of freedom, the JT instabilities do not contribute to the formation of the noncollinear 3 ordered magnetic state, which is present also in the undistorted cubic phase.
V Summary and Conclusions
In conclusion, in this study we have addressed the issue of parametrization of the LSDA+ approximation for noncollinear spin systems and explained the origin of the canted 3 state in UO by combining an array of advanced computational techniques including the constrained random phase approximation to compute the Coulomb repulsion parameter and the Hund’s coupling parameter , thus rendering the LSDA++SOC fully ab initio, magnetic constrains to study the dependence of the total energy on the direction of the spin moment, adiabatic propagation of the occupation matrix to avoid the multiple minima problem in constructing the accurate magnetic energy landscape, and two different effective Hamiltonians to extract from the ab initio data the value of the spinorbit interaction (achieved using an atomiclike relativistic Hamiltonian for orbitals) and the quadrupolequadrupole exchange interactions (using an effective pseudospin Hamiltonian).
The outcome of our study is threefold. First, we have derived an invariant orbital and spindependent formalism for the LSDA+ model for noncollinear spins that involves a spin and orbital convolution of the density matrix, and have shown that, unlike in collinear LSDA+, the noncollinear LSDA+ potential and double counting correction depend on one parameter only, and is independent on Hund’s coupling parameter . Second, our data suggest that the spinorbit interaction parameter in UO is as large as 1.45 eV, which is among the largest SOC strengths ever observed in a bulk material, and this represents the origin of many exotic physical phenomena emerging from the unusually intricate interplay between the spin, charge and orbital degrees of freedom, explicated by the formation of a multipolar magnetic state associated with tilted orbital ordering in the orbital manifold. Finally, we have uncovered the role of dipoledipole and quadrupolequadrupole spin interactions in the formation of the noncollinear 3 state and ruled out JahnTeller distortions as a factor in stabilizing the 3 magnetic ordering. The most relevant energy scales defining the properties of UO are summarized in Table 7.