Hubbard-like Hamiltonians for interacting electrons in s, p and d orbitals

Hubbard-like Hamiltonians are widely used to describe on-site Coulomb interactions in magnetic and strongly-correlated solids, but there is much confusion in the literature about the form these Hamiltonians should take for shells of p and d orbitals. This paper derives the most general s, p and d orbital Hubbard-like Hamiltonians consistent with the relevant symmetries, and presents them in ways convenient for practical calculations. We use the full configuration interaction method to study p and d orbital dimers and compare results obtained using the correct Hamiltonian and the collinear and vector Stoner Hamiltonians. The Stoner Hamiltonians can fail to describe properly the nature of the ground state, the time evolution of excited states, and the electronic heat capacity.


I. INTRODUCTION
The starting point for any electronic structure calculation is the Hamiltonian that describes the dynamics of the electrons. Correlated electron problems are often formulated using Hamiltonians that involve only a small number of localized atomic-like orbitals and use an effective Coulomb interaction between electrons that is assumed to be short ranged, and thus confined to within atomic sites. However, there is widespread variation in the literature over the form chosen for these Hamiltonians, and they often do not retain the correct symmetry properties. Here we derive the most general Hamiltonians of this type that are rotationally invariant and describe electrons populating s, p and d orbitals.
Despite often missing terms, multi-orbital Hamiltonians with on-site Coulomb interactions have been used successfully to describe strong correlation effects in systems with narrow bands near the Fermi energy. Applications include studies of metal-insulator transitions 1 , colossal magnetoresistant manganites 2 , and d band superconductors such as copper oxides 3 and iron pnictides 4 . The Hubbard Hamiltonian has also been used to simulate graphene 5,6 , where both theory and experiment have found magnetic ordering in nanoscale structures 7,8 . However, the lack of consistency in the forms used for the Hamiltonians is problematic. Indeed, according to Dagotto 2 , "the discussions in the current literature regarding this issue are somewhat confusing". The Hubbard Hamiltonian 9 is both the simplest and best known Hamiltonian of the form we are considering, and is valid for electrons in s orbitals. Hubbard was not the first to arrive at this form; two years earlier, in 1961, Anderson wrote down a very similar Hamiltonian 10 for electrons in s orbitals with on-site Coulomb interactions. Furthermore, in the appendix of the same paper, the Hamiltonian was extended to two orbitals, introducing the on-site exchange interaction, J. Unfortunately, this extra term is not rotationally invariant in spin space. In 1966, Roth 11 also considered two orbitals and two parameters, producing a Hamiltonian that is rotationally invariant in spin space. In 1969 Caroli et al. 12 corrected Anderson's Hamiltonian for two orbitals, again making it rotationally invariant in spin space. However, the Hamiltonians of Roth and Caroli et al. do not satisfy rotational invariance in orbital space. This was corrected by Dworin and Narath in 1970 13 , who produced a Hamiltonian that is rotationally invariant in both spin and orbital space. They generalized the multi-orbital Hamiltonian to include all 5 d orbitals, although they only included two parameters: the on-site Hartree integral, U , and the onsite exchange integral, J. The next contribution to the multi-orbital Hamiltonian was by Lyon-Caen and Cyrot 14 , who in 1975 considered a two-orbital Hamiltonian for the e g d orbitals. They introduced U 0 − U = 2J, where U 0 = V αα,αα is the selfinteraction term, U = V αβ,αβ is the on-site Hartree integral, and J = V αβ,βα is the on-site exchange integral, important for linking these parameters. Here V is the Coulomb interaction between electrons, and α and β are indices for atomic orbitals. Lyon-Caen and Cyrot also included the pair excitation, Jĉ † I1,↑ĉ † I1,↓ĉ I2,↓ĉI2,↑ , wherê c † andĉ are fermionic creation and annihilation operators. The new term moves a pair of electrons from a completely occupied orbital, 2, on site I, to a completely unoccupied orbital, 1, on site I. Such double excitations are common in the literature [15][16][17][18] . In 1978, Castellani et al. 19 wrote down the three-orbital Hamiltonian for the t 2g d orbitals 20 in a very clear and concise way: its form is exactly what we find here for p orbital symmetry. They too make note of the equation U 0 − U = 2J and include the double excitation terms. The next major contribution was by Oleś and Stollhoff 21 , who introduced a d orbital Hamiltonian with three independent parameters. The third parameter, ∆J, represents the difference between the t 2g and the e g d orbital exchange interactions; they estimated ∆J to be 0.15J. In our discussion of d orbitals we adopt their notation. Unfortunately, their Hamiltonian is not rotationally invariant in orbital space. In 1989, Nolting et al. 22 referenced the Hamiltonian proposed by Oleś and Stollhoff but discarded the double excitation terms and ignored the ∆J term. Their Hamiltonian will be called the vector Stoner Hamiltonian in this paper. Since then, most papers have used a Hamiltonian similar to that of Nolting et al., with only two parameters; see for example references 23-25. The vector Stoner Hamiltonian is often simplified by replacing the rotationally invariant moment-squared operator,m 2 =m ·m, bŷ m 2 z , yielding the collinear Stoner Hamiltonian. The multi-orbital Hamiltonians derived here consist of one-particle hopping (inter-site) integrals and twoelectron on-site Coulomb interactions.
The on-site Coulomb interactions for s, p and d orbital symmetries are presented in a clear tensor notation, as well as in terms of physically meaningful, rotationally invariant (in both orbital and spin space), operators that include all of the terms from the previous papers as well as additional terms of the same order of magnitude that have not been previously considered. The d orbital Hamiltonian presented in section II D corrects the d orbital Hamiltonian proposed by Oles and Stollhoff 21 by restoring full rotational invariance in orbital space. For an s orbital Hamiltonian, the on-site Coulomb interactions may be described using one independent parameter. Two independent parameters are required for p orbital symmetry, and three for d orbital symmetry. Furthermore, the method presented can be extended to f and g orbital symmetry and generalized to atoms with valence orbitals of multiple different angular momenta. The use of a restricted basis set and the neglect of the two-electron inter-site Coulomb interactions means that screening effects are missing. Thus, even though the parameters that define the s, p and d Hamiltonians appear as bare on-site Coulomb integrals, screened Coulomb integrals should be used when modelling real systems; this significantly reduces the values of the parameters. Although it is possible to find the on-site Coulomb integrals using tables of Slater-Condon parameters 26,27 or the closed forms for the integrals of products of three spherical harmonics found by Gaunt 28 and Racah 29 , the point of this paper is to remove that layer of obscurity and present model Hamiltonians written in a succinct form with clear physical meaning. The Slater-Condon parameters and Gaunt integrals were not used in the derivations of the on-site Coulomb interaction Hamiltonians presented here; however, they did provide an independent check to confirm that the tensorial forms derived here are correct. The link between the spherical harmonics Coulomb integrals, the Slater-Condon parameters and the Racah parameters is explained clearly in 30 . The transformation to cubic harmonics is straightforward; an example of which has been included for the p-orbitals in Section II C. We note that Rudzikas 31 has derived rotationally invariant Hamiltonians for p and d orbital atoms. However, his derivation is different from ours and his d orbital Hamiltonian is partially expressed in terms of tensors; the form we present in this paper is written in terms of physically meaningful operators. The p and d symmetry Hamiltonians are described in Section II (see Appendix B for full derivations) and compared with the vector Stoner Hamiltonian. In Section III we show the ground state results for both our Hamiltonian and the vector Stoner Hamiltonian for p and d orbital dimers. We find that, although the results are often similar, there are important qualitative and quantitative differences. In Section IV we show the differences that emerge when evolving a starting state for a p orbital dimer under our Hamiltonian and the vector Stoner Hamiltonian. We also calculate a thermodynamic property, the electronic heat capacity of a p orbital dimer, and obtain significantly different results when using our Hamiltonian and the vector Stoner Hamiltonian. The eigenstates of all Hamiltonians are found using the full configuration interaction (FCI) method as implemented in the HANDE code 32 .

A. General form of the Hamiltonian
We start with the general many-body Hamiltonian of electrons in a solid, expressed in second-quantized form using a basis of localized one-electron orbitals and ignoring terms that do not include electronic degrees of (1) Hereĉ † Iα,σ andĉ Iα,σ are creation and annihilation operators, respectively, for an electron in orbital α on site I with spin σ. Upper case Roman letters such as I, J, K and L refer to atomic sites, the lower case Greek letters α, β, χ and γ represent localized orbitals, and σ and ζ indicate spin (either up or down). The one-electron contributions to the Hamiltonian are encapsulated in the t IαJβ term that includes the electronic kinetic energy and electron-nuclear interaction. The electron-electron interaction terms are represented by the Coulomb matrix elements V IαJβ,KχLγ . We now make two approximations: we only retain electron-electron Coulomb interactions on-site and we restrict the basis to a minimal set of localized angular momentum orbitals. Thus we have just one s orbital per site for the s-symmetry Hamiltonian, three p orbitals per site for the p-symmetry Hamiltonian, and five d orbitals per site for the d-symmetry Hamiltonian. The Hamiltonian can then be written as follows: where V I αβ,χγ = V IαIβ,IχIγ . To improve readability, from this point onwards we drop the site index, I, from the Coulomb integrals and creation and annihilation operators as we shall always be referring to on-site Coulomb interactions. The on-site Coulomb interaction, V αβ,χγ , is a rotationally invariant tensor if the atomic-like orbitals are angular momentum eigenstates. To demonstrate this we note that applying the rotation operatorR to an atomic-like orbital φ α of angular momentum giveŝ where d βα (R) is the (2 + 1) × (2 + 1) matrix that corresponds to the rotationR in the irreducible representation of angular momentum . By making the change of variable x =Rr and x =Rr , it is straightforward to show that which demonstrates that V αβ,χγ is indeed a rotationally invariant tensor.
B. The one band Hubbard model: s orbital symmetry The Hubbard Hamiltonian is applicable only to one band models where the one orbital per atom has s symmetry. This gives a simple form for the on-site Coulomb interaction tensor, U 0 = V αα,αα ; as there is only one type of orbital, there is only one matrix element. We can use this result to simplify the Coulomb interaction part of Eq. (2) and express it in terms of the electron number operator or the magnetic moment vector operator, defined as follows: where σ are the Pauli spin matrices and the sum over α has only one term as there is only one spatial orbital for the s case. Hence or, equivalently, where we have use the normal ordering operator, ::, to remove self interactions. The action of this operator is to rearrange the creation and annihilation operators such that all the creation operators are on the left, without adding the anticommutator terms that would be required to leave the product of operators unaltered; if the rearrangement requires an odd number of flips, the normal ordering also introduces a sign change. For example :n 2 : = αβσζ :ĉ † α,ζĉα,ζĉ † β,σĉβ,σ :, If the normal ordering operator is not used in Eqs. (7) and (8) then additional one electron terms have to be subtracted to remove the self interaction. The mean-field form of this Hamiltonian has been included in Appendix A 1.

C. The multi-orbital model Hamiltonian: p orbital symmetry
Consider the case where the local orbitals, α, β, χ and γ appearing in Eq. (4) are real cubic harmonic p orbitals with angular dependence x/r, y/r, and z/r. In this case the rotation matrices, d α α (R), are simply 3 × 3 Cartesian rotation matrices. This makes V αβ,χγ a rotationally invariant fourth-rank Cartesian tensor, the general form of which is 34 , where U = V αβ,αβ , J = V αβ,βα and J = V αα,ββ , with α = β. By examination of the integral in Eq. (4), we note that the Coulomb tensor has an additional symmetry when the orbitals are real (as are the cubic harmonic p orbitals), namely V αβ,χγ = V χβ,αγ and V αβ,χγ = V αγ,χβ , and hence J must be equal to J . Thus we find This recovers the well known equation 14,19 , and shows that the most general cubic p orbital interaction Hamiltonian is defined by exactly two independent parameters 35 . In Appendix B 1 we give a full derivation of Eq. (11) using representation theory. In passing we observe that symmetric fourth-rank isotropic tensors can be found in other areas of physics, such as the stiffness tensor in isotropic elasticity theory 36 where λ is the Lamé coefficient and µ is the shear modulus. Transforming Eq. (11) into a basis of complex spherical harmonic p orbitals with angular dependence of the form and has the following values The result of applying this transformation is where m, m , m and m are spherical harmonics indices going from −1 to +1, U = V 10,10 , and J = V 10,01 . Using Eq. (11) as a convenient starting point, it is straightforward to rewrite the on-site Coulomb Hamiltonian in terms of rotationally invariant operators: 31 where the vector angular momentum operator for p orbitals on site I is defined to bê and µβα is the three-dimensional Levi-Civita symbol. An equivalent expression is where the final term corresponds to on-site electron Eq. (16) embodies Hund's rules for an atom. By making the substitutionm = 2Ŝ, we see that the spin is maximized first (prefactor −2J) and then the angular momentum is maximized (prefactor − 1 2 J). The mean-field form of the above p orbital Hamiltonian is given in Appendix A 2.

D. The multi-orbital model Hamiltonian: d orbital symmetry
The on-site Coulomb interaction for cubic harmonic d orbitals can be expressed as where ξ is a five-component vector of traceless symmetric 3 × 3 transformation matrices defined as follows We use the convention that index numbers (1,2,3,4,5) correspond to the d orbitals (3z 2 −r 2 , zx, yz, xy, x 2 −y 2 ). See Appendix B 2 for a derivation of Eq. (21) using representation theory. We have chosen U to be the Hartree term between the t 2g orbitals, J to be the average of the exchange integral between the e g and t 2g orbitals, and ∆J to be the difference between the exchange integrals for e g and t 2g . That is, This choice of parameters is the same as was used by Oleś and Stollhoff 21 ; Eq. (20) generalizes their result and restores rotational invariance in orbital space. Note that, in the mean field, the on-site block of the density matrix in a cubic solid is diagonal; this will cause the mean-field version of Eq. 22) to reduce to that of Oleś and Stollhoff 21 . However the presence of interfaces, vacancies or interstitials will break the cubic symmetry make the on-site blocks of the density matrix non-diagonal. Therefore the use of the complete Hamiltonian is recommended in mean-field calculations for systems that do not have cubic symmetry. Rewriting Eq. (20) in terms of rotationally invariant operators giveŝ whereQ 2 is the on-site quadrupole operator squared. The quadrupole operator for a single electron is defined asQ The quadrupole operator for a system of many electrons is obtained by summing the operators for each electron.
Since each term in the sum acts on the coordinates of one electron only, the operatorsL µ andL ν appearing in that term measure the angular momentum of one electron only, not the whole system. More detail on theQ 2 operator is included in Appendix C. Eq. (22) is similar in form to that found in Rudzikas 31 , the difference being that our Hamiltonian is expressed in terms of physically meaningful operators. The mean-field form of this Hamiltonian is given in Appendix A 3.

E. Comparison with the Stoner Hamiltonian
The Stoner Hamiltonian for p and d orbital atoms, as generally defined in the literature, has the following form for the on-site Coulomb interaction: 25 where J, in this many-body context, is the Stoner parameter, I. However, in the mean field, the Stoner parameter should not be taken to be equal to J due to the self-interaction error; see reference 38 . This form clearly breaks rotational symmetry in spin space, which can be restored by substituting :m 2 z : with :m 2 : In this paper we will compare our Hamiltonians, Eqs. (18) and (22), with Eq. (25), which we shall call the vector orm 2 Stoner Hamiltonian. By starting with the vector Stoner Hamiltonian and working backwards, it is possible to find the corresponding tensorial form of the Coulomb interaction: This is very similar to the tensorial form for the p-case Coulomb interaction, Eq. (11), except that it is missing a term, Jδ αβ δ γχ . The lack of this term means that the vector Stoner Hamiltonian does not respect the symmetry between pairs αχ and βγ evident from the form of the integral in Eq. (4), i.e., Vm 2 Stoner χβ,αγ = Vm 2 Stoner αβ,χγ . As we shall see later, this omission changes the computed results.

A. Simulation setup
Here we investigate the eigenstates of dimers using both our Hamiltonians, Eqs. (18) and (22), and the vector Stoner Hamiltonian, Eq. (25). We employ a restricted basis with a single set of s, p and d orbitals for the s, p and d Hamiltonians respectively. We perform FCI calculations using the HANDE 32 computer program to find exact eigenstates of the model Hamiltonians. The Lanczos algorithm was used to calculate the lowest 40 wavefunctions for our Hamiltonian, more than sufficient to find the correct ground state wavefunction, while the full spectrum was calculated for the vector Stoner Hamiltonian. The single particle contribution to the Hamiltonian for all models (the hopping matrix) is defined for a dimer aligned along the z-axis: in this case it is only non-zero between orbitals of the same type, whether on the same site or on different sites. The sigma hopping integral, |t σ |, was set to 1 to define the energy scale. The other hopping matrix elements follow the canonical relations suggested by Andersen 39 , and we use the values of Paxton and Finnis 40 : ppσ : ppπ = 2 : −1 for p orbitals, and ddσ : ddπ : ddδ = −6 : 4 : −1 for d orbitals.
The directionally averaged magnetic correlation between the two sites has been calculated for the ground state. If positive it shows that the atoms are ferromagnetically correlated (spins parallel on both sites) and if it is negative it shows that the atoms are antiferromagnetically correlated (spins antiparallel between sites). The correlation between components of spin projected onto a direction described by angles θ and φ, in spherical coordinates, is defined as follows where the :: is the normal ordering operator, used to remove self-interaction, and σ θφ ξξ = σ z ξξ cos θ + σ x ξξ sin θ cos φ + σ y ξξ sin θ sin φ. The magnetic correlation averaged over solid angle is B. Ground state: p orbital dimer The ground states of our Hamiltonian and the vector Stoner Hamiltonian have been found to be rather similar for 2, 6 and 10 electrons split over the p orbital dimer, but qualitatively different for 4 and 8 electrons. The magnetic correlation between the two atoms for 4 electrons is shown in Fig. 1. The symmetries of the wavefunctions 41 are indicated on the graph by the notation 2S+1 Λ ± u/g : the ± is only used for L z = 0 (i.e. Σ states) and corresponds to the sign change after a reflection in a plane parallel to the axis of the dimer; the u/g term refers to ungerade (odd) and gerade (even) and corresponds to a reflection through the midpoint of the dimer; Λ is the symbol corresponding to the total value for L z (e.g. L z = 1 is Π, L z = 2 is ∆, L z = 3 is Φ, L z = 4 is Γ); and S is the The bottom graph has a region with symmetry 3 Σ − g extending a long way up the J axis, which is not present in the ground state of our Hamiltonian. The bottom graph also includes a region with two degenerate states with symmetries 1 ∆ g and 1 Σ + g ; this degeneracy is broken when our Hamiltonian is used.
total spin. The differences between the ground states of our Hamiltonian and the vector Stoner Hamiltonian shown in Fig. 1 are as follows: the ground-state wavefunction of the vector Stoner Hamiltonian has a region with 3 Σ − g symmetry extending far up the J axis, which is not present for our Hamiltonian; the region with 1 ∆ g symmetry is doubly degenerate for our Hamiltonian, as total L z = ±2, whereas for the Stoner Hamiltonian it is triply degenerate, as it is also degenerate with the state of symmetry 1 Σ + g (this state appears at a higher energy for our Hamiltonian).

C. Ground state: d orbital dimer
The ground states of our Hamiltonian and vector Stoner Hamiltonian have been found to be rather similar for 2, 4, 6, 14, and 18 electrons split over the d orbital dimer when ∆J is small. The simulations for 10 electrons were not carried out as they are too computationally expensive. Qualitative differences were found for 8 and 12 electrons; for an example see Fig. 2, which shows the magnetic correlation between two atoms in the ground state of the d-shell dimer with 6 electrons per atom. From Fig. 2 we see that the results for our Hamiltonian with ∆J/|t| = 0.0 (top graph) and with a small value of ∆J/|t| = 0.1 (middle graph) are qualitatively different from those for the vector Stoner Hamiltonian (bottom graph). The vector Stoner Hamiltonian makes the 1 Γ g and 1 Σ + g states degenerate, which is not the case for our Hamiltonian. The largest differences are between the ∆J/|t| = 0.1 graph and the vector Stoner Hamiltonian: new regions with symmetry 1 Σ + g , 3 Φ g and 5 Σ − g appear, and the region with symmetry 3 Σ − g almost disappears. This shows that the inclusion of the quadrupole term in Eq. (22) can make a qualitative difference to the ground state.

IV. EXCITED STATE RESULTS
Differences between our Hamiltonian and the vector Stoner Hamiltonian are also observed for excited states.
Here we present three examples, one demonstrating explicitly the effect of the pair hopping, one showing more general differences in the excited states through a calculation of the electronic heat capacity as a function of temperature, and one showing a difference in the spin dynamics of a collinear and a non-collinear Hamiltonian. The full spectrum of eigenstates was calculated for all of these examples.

A. Pair hopping
The pair hopping term, J αβσσ ĉ † α,σĉ † α,σ ĉ β,σ ĉ β,σ = J αβ : (n αβ ) 2 :, is found in both our p and d orbital Hamiltonians, but is absent from the vector Stoner Hamiltonian. To demonstrate the effect of the pair hopping term we initialize the wavefunction in a state with two electrons in the x orbital on site 1 of a p-atom dimer. Figure 3 shows the evolution of the wavefunction in time using our Hamiltonian and the vector Stoner Hamiltonian, with U = 5.0|t| and J = 0.7|t|. For our Hamiltonian the pair of electrons in orbital x on atom 1 hops to the y and z orbitals on atom 1 (mediated by the pair hopping term) as well as between the two atoms (mediated  , the x orbital on atom 1, of a p-atom dimer with U = 5.0|t| and J = 0.7|t|. We see that the two electrons in 1x are able to hop into 1y and 1z when the wavefunction is evolved using our Hamiltonian but not when it is evolved using the vector Stoner Hamiltonian. by the single electron hopping term). The result is that there is a finite probability of finding the pair of electrons in any of the x, y and z orbitals on either atom. However, for the vector Stoner Hamiltonian, the two electrons in the x orbital on site 1 are unable to hop into the y and z orbitals on atom 1; they are only able to hop to the x orbital on atom 2. This means tha there is no possibility of observing the pair of electrons in the y or z orbitals on either atom.

B. Heat capacity
The heat capacity is calculated from where T is the temperature, k B is Boltzmann's constant, ε i is a many-electron energy eigenvalue, and N a is the number of atoms. The result for both Hamiltonians as a function of temperature for four electrons split over two p-shell atoms with U = 5.0|t| and J = 0.7|t| is shown in Figure 4. There is a qualitative difference between the results for our Hamiltonian (the solid line) and the vector Stoner Hamiltonian (dashed line). The energies of the eigenstates of our Hamiltonian are more spread out than those of the Stoner Hamiltonian. This is due to the inclusion of the pair hopping term, which causes states with on-site paired electrons to rise in energy (by ∼ J/|t|). The unphysical reduction of the heat capacity to zero as the temperature tends to infinity is a consequence of the use of a restricted basis set.

C. Spin dynamics
Here we show how the collinear Stoner Hamiltonian can give rise to unphysical spin dynamics. We consider an S z = +1 triplet with both electrons on one of the two p orbital atoms, |Ψ = |T +1 . Written as linear combinations of two-electron Slater determinants, the three states in the triplet are |Ψ rot is still an eigenstate, with eigenvalue U −J, of both the vector Stoner Hamiltonian and our Hamiltonian, but it is no longer an eigenstate of the collinear Stoner Hamiltonian. The wavefunction |Ψ rot evolves in time as

and our Hamiltonian
, where τ is time. If we take the expectation value of the magnetic moment on each site using the collinear Stoner Hamiltonian, we find In contrast, using the vector Stoner Hamiltonian or our p-shell Hamiltonian, the expectation value of the magnetic moment is independent of time, equal to (1, 0, 0) in vector format. This demonstrates that calculations of spin dynamics using the collinear Stoner Hamiltonian can give rise to unphysical oscillations of the magnetic moments.

V. CONCLUSION
We have established the correct form of the multi-orbital model Hamiltonian with on-site Coulomb interactions for atoms with valence shells of s, p and d orbitals. The methodology used may be extended to atoms with f and g shells, and to atoms with valence orbitals of several different angular momenta. The results presented show that there are important differences between our p-and d-shell Hamiltonians and the vector Stoner Hamiltonian. The vector Stoner Hamiltonian misses both the pair hopping term, which is present in our p and d orbital Hamiltonians, and the quadrupole term, which is present in our d orbital Hamiltonian. The pair hopping term pushes states with pairs of electrons up in energy, whereas the magnetism term pulls states with local magnetic moments down in energy. The pair hopping term has the largest effect on the ground state for p orbitals with 2 and 4 electrons per atom and for d orbitals with 4 and 6 electrons per atom, for J/|t| < 2. This is because the number of possible determinants with paired electrons on site is large for these filling factors and the low lying states can be separated based on this. At values of J/|t| > 2 the magnetism becomes the dominant effect upon the selection of the ground state and the difference between the ground state of our Hamiltonian and that of the vector Stoner Hamiltonian becomes small. However, the differences are rather more pronounced for the excited states. This is evidenced by the hopping of pairs of electrons between orbitals on the same site, which is allowed by our Hamiltonian but not by the vector Stoner Hamiltonian, and in differences in the electronic heat capacity as a function of temperature. We also find clear evidence that the collinear Stoner Hamiltonian is inappropriate for use in describing spin dynamics as it breaks rotational symmetry in spin space. We would expect similar problems when using collinear time dependent DFT simulations to model spin dynamics. Application of the mean-field approximation 42 to the onsite Coulomb interaction part of the Hubbard Hamiltonian with s orbital symmetry, Eq. (7), yields the following,V for which the total energy is where the double counting correction has been included. Equivalently, applying the mean-field approximation to Eq. (8) yields the following, for which the total energy is where again the double counting correction has been included.

The multi-orbital model Hamiltonian: p orbital symmetry
Application of the mean-field approximation 42 to the model Hamiltonian with p orbital symmetry, Eq. (18), yields the following, for which the total energy is where the double counting correction has been included.

The multi-orbital model Hamiltonian: d-orbital symmetry
Application of the mean-field approximation 42 to the model Hamiltonian with d-orbital symmetry, Eq. (22), yields the following, αβ n αβ n βα + n αβ n αβ for which the total energy is where the double counting correction has been included. As a precursor to our treatment of the d case in Appendix B 2, we use representation theory 43 to derive Eq. (11) for cubic harmonic p orbitals. First we define the following irreducible objects: "0" is a scalar, "1" is a threedimensional vector, "2" is a 3 × 3 symmetric traceless second-rank tensor, "3" is a third-rank tensor consisting of three "2"s, and "4" is a fourth-rank tensor consisting of three "3"s. The Coulomb interaction V αβ,γδ transforms as a tensor product of four "1"s, 1 ⊗ 1 ⊗ 1 ⊗ 1, but we know also that it must be isotropic and that from the above irreducible objects only the "0" is isotropic. Objects such as 1 ⊗ 1, represented as A ij , may be expanded just as in the addition of angular momentum: where "0" is a scalar, s = ij A ij δ ij , "1" is a vector, v i = jk A jk ijk , and "2" is a traceless symmetric ten- The reader may be more familiar with combinations of angular momentum. In the language of angular momentum, the transformation properties of a tensor product of two objects of angular momentum with l = 1 are described by a 9 × 9 matrix, which is block diagonalisable into a 1 × 1 matrix (which describes the transformations of an l = 0 object), a 3 × 3 matrix (which describes the transformations of an l = 1 object), and a 5 × 5 matrix (which describes the transformations of an l = 2 object). In the angular momentum representation the matrix elements are complex. Here we are using cubic harmonics and the matrix elements are real. Similarly 1 ⊗ 1 ⊗ 1 ⊗ 1 may be expanded out as: The only isotropic, "0", objects arise from 0⊗0, 1⊗1 and 2 ⊗ 2. A general isotropic three-dimensional fourth-rank tensor, T ijkl , therefore has three independent scalars that are found from symmetric and antisymmetric contractions of T ijkl . We are interested in the form of the p orbital on-site Coulomb interactions, V αβ,χγ , which have an additional symmetry, such that they remain unchanged under exchange of α and χ and under exchange of β and γ. We therefore only require the symmetric scalars that arise from the 0 ⊗ 0 and 2 ⊗ 2 contractions to describe V αβ,χγ : V αβ,χγ = s 0 δ αχ δ βγ + s 2 (δ αγ δ βχ + δ αβ δ χγ ) , which is symmetric under exchange of α and χ or β and γ; where This is equivalent to Eq. (11), where s 0 = U and s 2 = J.

d orbital symmetry
To find the isotropic five-dimensional fourth-rank tensor that describes the on-site Coulomb interactions for d orbital symmetry, we find it convenient to map from a fivedimensional fourth-rank tensor to a three-dimensional eighth-rank tensor. We do this by replacing the fivedimensional vector of d orbitals by the three-dimensional traceless symmetric B matrix of the d orbitals: (B5) We refer to the elements of B as B ab , and use the subscripts a, b, c, d, s, t, u and v as the indices for irreducibles of rank 2 or higher. The transformation between the irreducible B ab and the cubic harmonic d orbitals is where ξ is a traceless, symmetric transformation matrix, defined in Eq. (21), π is a normalisation factor, R d is the radial function for a cubic d orbital, and S is the spin function. The orbital indices α, β, χ and γ run over the five independent d orbitals, whereas the irreducible indices a, b, c, d, s, t, u and v run over the three Cartesian directions. The mapping between the five-dimensional fourth-rank tensor and the three-dimensional eighth-rank tensor is as follows where this equation defines the isotropic threedimensional eighth-rank tensor for the on-site Coulomb integrals, V ab,cd,st,uv , and we have dropped the complex conjugates on the ξ matrices as they are real. The B matrix is a three-dimensional traceless symmetric "2" and thus contains 5 independent terms. As we show below, it is possible to represent the isotropic threedimensional eighth-rank tensor as a list of quadruple Kronecker deltas, with 5 independent parameters in general and 3 independent parameters for the on-site Coulomb interactions. The isotropic three-dimensional eighthrank tensor is then converted back into an isotropic fivedimensional fourth-rank tensor. For an independent verification of the form of the isotropic three-dimensional eighth-rank tensor see Ref. 44. The power of representation theory is that one can follow an analogous procedure for shells of higher angular momentum, using a "3" to represent the seven f orbitals, a "4" to represent the nine g orbitals, and so on. One can also construct interaction Hamiltonians for atoms with important valence orbitals of several different angular momenta.
We now return to the case of a d shell and explain the procedure used to find the general isotropic Coulomb Hamiltonian in more detail. We start by representing the isotropic three-dimensional eighth-rank tensor, V ab,cd,st,uv , as a 2 ⊗ 2 ⊗ 2 ⊗ 2. Proceeding as we did for the p case, we write The terms in Eq. (B8) that can generate a rank 0 (a scalar) are 0 ⊗ 0, 1 ⊗ 1, 2 ⊗ 2, 3 ⊗ 3 and 4 ⊗ 4. This implies that a general isotropic three-dimensional eighthrank tensor is defined by five independent parameters. We know that, for V αβ,χγ , there exists a symmetry between the pairs αχ and βγ. It follows that V ab,cd,st,uv has a symmetry between the pairs (ab)(st) and (cd)(uv). We therefore require only the even contractions, 0 ⊗ 0, 2 ⊗ 2 and 4 ⊗ 4, reducing the five parameters to three 45 ; this is the same number proposed in Chapter 4 of Ref. 35.
To get the scalars from V ab,cd,st,uv we have to contract it. There are six possible contractions of V ab,cd,st,uv , outlined in Table I, of which only five are independent. Note that one does not contract within the indices ab, cd, st or uv because B ab is traceless, hence ab δ ab B ab = 0. By examining the diagrammatic contractions from Table I We find that the last three of the transformations of the contractions, although rotationally invariant in orbital space, are not expressible in terms of Kronecker deltas of α, β, χ and γ. This is to be expected as there exist terms in V αβ,χγ that are non-zero and have more than two orbital indices, e.g. V (3z 2 −r 2 )(xy),(yz)(xz) or V (3z 2 −r 2 )(xy),(xy)(x 2 −y 2 ) . We can now write down V αβ,χγ concisely, where we have defined U to be the Hartree term between the t 2g orbitals on-site, U = V (zx)(yz),(zx)(yz) , J to be the average of the exchange integral between the e g and t 2g orbitals on-site, J = 1 2 V (zx)(yz),(yz)(zx) + V (3z 2 −r 2 )(x 2 −y 2 ),(x 2 −y 2 )(3z 2 −r 2 ) , and ∆J to be the difference between the exchange integrals for e g and t 2g , ∆J = V (3z 2 −r 2 )(x 2 −y 2 ),(x 2 −y 2 )(3z 2 −r 2 ) − V (zx)(yz),(yz)(zx) . We could equally well have chosen to use the sum of the fourth and sixth contractions instead of the fifth contraction, although with different coefficients.
Appendix C: Angular momentum and quadrupole operators

Angular momentum
The angular momentum operators are most naturally expressed in spherical polar coordinates, even when using them to operate on cubic harmonics. Their forms are: where we have seth to 1 for simplicity. When applying the above operators to the on-site cubic harmonic p orbitals (we have dropped the site index for clarity) it is straightforward to show that A general one particle operatorÔ may be expressed in terms of creation and annihilation operators as follows: Substituting Eq. (C4) into Eq. (C5) yieldŝ where µ is a Cartesian direction and α and β are cubic harmonic p orbitals (x, y or z). The d case is slightly more complicated, but is greatly simplified by the use of Eq. (B6). Applying the angular momentum operators to the B matrix reveals that witĥ L µ |φ ασ = iN d jks S(σ) R d (r) r 2 ξ αjk ( µjs |B sk + µks |B js ) .

(C8)
By finding the expectation value and substituting into Eq. (C5) we obtain