Ligand molecular orbitals (MOs) have been recognized as important players in the physics of transition-metal compounds since the introduction of the Zhang-Rice singlet back in 1988 . They have recently gained a renewed interest after the concept of oxygen MOs strongly coupled to lattice degrees of freedom in a polaronic way was used to describe rare-earth nickelates   and superconducting bismuthates   . Ligand MOs are particularly important in negative charge-transfer and hole-doped charge-transfer insulators, where often they are the orbitals that end up being occupied by the (self-)doped holes and therefore have a direct impact on the system’s low-energy properties . Among this important class of materials are such vigorously discussed but still controversial systems as superconducting cuprates   and bismuthates   , rare-earth nickelates , and transition-metal oxide based Li-ion battery cathode materials  . Since this topic is expected to generate only more interest in the future, the current paper will present a band theory perspective on some important aspects of MOs in complex oxides, including the mathematical construction of MO-based basis sets (Section 2), the non-orthogonality problem (Section 3), their usage in effective models in the presence of strong electron-phonon coupling (Section 4), as well as the problem of a proper treatment of strong local correlations on MOs (Section 5); conclusions will be given in Section 6. Each part of the discussion will include a detailed example of a real material.
2. Projecting Electronic States onto Molecular Orbitals (NaNiO2 as an Example)
Let us consider the prototypical negative charge-transfer insulator NaNiO2 as our first example. Its monoclinic C2/m crystal structure consists of triangular layers of edge-sharing cooperatively elongated NiO6 octahedra intercalated with Na+ ions, with one formula unit per primitive unit cell . Although here Ni has a formal valency of 3+, spectroscopic studies clearly show the abundant presence of oxygen holes which suggest an actual electronic configuration of Ni2+L, where L is an oxygen hole . This three-hole state has a net spin of ½ and would formally correspond to Ni3+ ( ) with spin 1/2, which is quite contrary to what one would expect from Hund’s rules if it really were Ni3+. Due to the strong hybridization between the Ni—eg orbitals and the oxygen MOs of respective symmetry, the oxygen holes select to occupy the O—(x2-y2) MOs, one hole per oxygen octahedron, and form a Zhang-Rice-like spin singlet state with the hole in the Ni—(x2-y2) orbital. This selective occupation of the (x2-y2) MOs by the holes, which results in the mentioned cooperative elongation of the NiO6 octahedra, constitutes a molecular orbital (or band) Jahn-Teller effect.
A way of verifying and visualizing this picture using band theory is to perform a projection of the NaNiO2 electronic structure on a MO-based basis set. We note, however, that band theory is unable to fully capture the multi-determinant character of the spin singlet state in question and will underestimate the energy
of the singlet as rather than , where Jpd is the spin exchange coupling
between interacting atomic d and molecular p states. Still, it provides a fair description of on-site symmetries and spin densities, worth exploring with the projection technique.
Figure 1(a) shows the calculated NaNiO2 density of states (DOS) projected onto the Ni—3d orbitals (left) and the oxygen molecular orbitals (right) for the majority and minority spin channels. This is a spin-polarized local density
Figure 1. Oxygen molecular orbitals in NaNiO2. (a) Projected densities of states (PDOS) per formula unit for the Ni—3d atomic orbitals and the oxygen MOs. (b) An isolated oxygen—(x2-y2) MO . The primitive cell comprising a single formula unit is marked with black lines. The short Ni—O bonds in the elongated octahedra are colored in a darker grey color than the long ones. (c) The majority spin band structure with the oxygen—(x2-y2) character highlighted in red. (d) An oxygen—(x2-y2) Bloch function at the point. In (a) and (c) the Fermi level is set to zero and marked with a dashed line. (e) Oxygen—pσ orbitals are in fact the , , and orbitals of the two inequivalent oxygens O1 and O2 in a unit cell, but each coming from an oxygen belonging to a different unit cell.
approximation + U (LDA + U)   calculation (on-site Coulomb repulsion U = 6 eV, Hund’s exchange interaction JH = 1 eV for Ni—3d electrons  ) performed, like the rest of the calculations presented in this paper, using the linearized augmented plane waves method implemented in the Wien2k package . A ferromagnetic order of Ni magnetic moments, both within and between the Ni planes, is assumed for simplicity. The majority spin hole has a mixed character of the Ni—(x2-y2) atomic orbital and the oxygen—(x2-y2) molecular orbital, which is a result of strong hybridization between these orbitals. This hole state is the anti-bonding combination of the two orbitals pushed up in energy above the Fermi level, while the bonding combination lands at about −5 eV below the Fermi level. An isolated oxygen—(x2-y2) molecular orbital centered at a Ni site R is shown in Figure 1(b), but the actual projection of the DOS [Figure 1(a)] or the k-resolved projection of the band structure [Figure 1(c)] are performed onto the Bloch functions constructed from :
where M indicates an MO character of the wave function and k is the crystal wave vector. To give a visual example, we show the oxygen—(x2-y2) MO Bloch function at the point [k = (0, 0, 0)] in Figure 1(d).
Generally speaking, MO Bloch functions can be obtained from the atomic orbital Bloch functions via a unitary transformation U(k):
Allowing for the k-dependence of the transformation matrix elements Uij(k) reflects the fact that an MO can be built from atomic orbitals nominally belonging to different unit cells. This is, for instance, the case for the oxygen—(x2-y2) MO shown in Figure 1(b), where the oxygen—pσ orbitals to the left of the Ni atom and those to its right belong to different unit cells. Figure 1(e) explains in further detail that the six oxygen pσ orbitals on a given NiO6 octahedron are in fact the , , and orbitals of the two inequivalent oxygens O1 and O2 in a unit cell, but each coming from an oxygen belonging to a different unit cell. In a way, using a k-dependent unitary transformation matrix U(k) broadens the technical definition of a unit cell as used in band-structure calculations for the purpose of setting up a Bloch basis set, where now atomic orbitals from the same atom can belong to different unit cells as in the case of NaNiO2.
Let us also note that if a MO character of an electronic state shows very little variation with k-vector this should be regarded as a sign of the robustness of the molecular nature of this state in real space. As one can see in Figure 1(c), this is the case for our example of NaNiO2, where indeed the majority spin hole state has a strong oxygen—(x2-y2) character at every k-vector.
Although in the charge- and negative charge-transfer systems it is very typical for the ligand holes to occupy MOs of the same symmetry as that of the cation atomic orbitals that they hybridize with [like the oxygen—(x2-y2) MO and the Ni—(x2-y2) atomic orbital in NaNiO2], there are a number of notable exceptions. For example, in the iron disulfide FeS2 (pyrite structure) and the recently discussed iron dioxide FeO2 , the ligand holes reside on the MOs formed on sulfur/oxygen dimers owing to strong intra-dimer hybridization between the sulfur/oxygen—pσ orbitals. The same is the case in the superoxides like KO2, where oxygen dimers have a net spin of 1/2 and which order magnetically below the ordering temperature.
As a final technical remark for this section, our MO projections are performed using atomic-like functions , which are the solutions of the Schrödinger equation within the muffin-tin sphere of atom α at the linearization energy and where are spherical harmonics, to construct molecular orbitals. This is in some contrast with the approach adopted in Wien2k, where atomic orbital projections are done onto the spherical harmonics inside muffin-tin spheres. In practice, the two approaches give very similar results for the projected DOS, but having also the radial part means that projections are being done on a true atomic-like orbital (or molecular orbital combination), which, in particular, enables their visualization in real space.
3. Non-Orthogonality (BaBiO3 as an Example)
For certain lattice structures—including those of the cuprates, with their characteristic two-dimensional square CuO2 planes, as well as of the bismuth and nickel perovskites—construction of MOs may also require their proper orthonormalization performed on top of the unitary basis transformation discussed above. In this regard, let us consider the example of the barium bismuth perovskite BaBiO3. In a perovskite lattice structure of the ABO3 type [Figure 2(a)], neighboring BO6 octahedra share their corners such that the O-B-O angle is equal or close to 180˚. As a result, the oxygen—pσ orbitals, typically involved in the formation of MOs, have to be shared between neighboring B cations. In the case of BaBiO3, the MO of most importance is oxygen—a1g, shown in Figure 2(a), due to its strong hybridization with the Bi—6s atomic orbital. However, because of the problem of a common oxygen—pσ orbital, neighboring oxygen—a1g MOs are non-orthogonal to each other. If no orthonormalization is done to such MOs, the Bloch functions constructed out of them will have a k-vector dependent overlap integral, , and may even completely vanish at certain k-points. As demonstrated pictorially in Figure 2(c) & Figure 2(d) and also shown in the LDA band structure calculation in Figure 2(b), this is the case for the oxygen—a1g MO in BaBiO3 at the point.
In order to be suitable for use in effective models, MOs of this kind need to be orthonormalized first. As was originally shown by Zhang and Rice , the k-vector dependent normalization factor ,
Figure 2. Oxygen molecular orbitals in BaBiO3. (a) An idealized (e.g., without octahedra’s rotations or bond disproportionation) cubic perovskite crystal structure of BaBiO3, with one formula unit per unit cell, featuring also an isolated oxygen—a1g MO . (b) Projection of the BaBiO3 electronic states onto oxygen MOs; the Fermi level is marked with a dashed line. (c) and (d) show oxygen—a1g MO Bloch functions at the [k = (0, 0, 0)] and R [k = (π, π, π)] points, respectively. Note that at the point the individual contributions completely cancel out due to their non-orthogonality, which results in vanishing oxygen—a1g MO weight at this k-vector in (b).
may be ill-defined at certain k-points. For the oxygen—a1g MO in BaBiO3, for example, it diverges at the point:
as obtained using the oxygen—a1g MO given by
where , i = x, y, z, are the Bloch functions of orthonormal oxygen—pσ orbitals and a is the cubic lattice constant. Let us note that an MO with contributions from both the oxygen—a1g MO and the Bi—6s atomic orbital ,
which was used by us and our co-authors in  for analyzing the BaBiO3 effective hopping parameters, does not have this divergence problem.
4. Comparison with Wannier Functions and Usage in Effective Models
It is important to recognize the fact that orthonormalization of a molecular orbital of any kind via multiplication of its corresponding Bloch function by converts this molecular orbital into a Wannier function . The same, of course, applies to atomic orbitals. It is also a common practice to adjust Wannier functions such that their corresponding eigenvalues match a selected set of bands obtained from an LDA(+U) calculation, which can be achieved through the following expansion in terms of self-consistently obtained Bloch eigenstates :
Here, can be either an atomic or a molecular orbital Bloch function and is a band index running from band n to band m. On top of this construction, one may also apply the maximal localization  and the disentanglement  procedures. This is finalized by orthonormalization [Equation (3)] to determine . The resulting Wannier functions are convenient to use in effective models since the basis set size can be kept minimal but at the same time the LDA(+U) bands would remain well reproduced by the model’s eigenstates. This approach is, for example, often used in regard to transition metal oxides in order to eliminate oxygen—p orbitals and derive a transition metal—d only based effective model. However, depending on how large and how strongly k-vector dependent is, such Wannier functions can be quite extended objects in real space with orbital contributions from many atomic shells. One disadvantage of using the Wannier functions [Equation (7)] in effective models is that most often it is not possible to fully control the degree of their spatial extension, even when attempts to optimize them through “maximal localization” have been made. This is a well-known problem for strongly localized atomic orbitals, like transition metal—3d, in Hubbard-type models where electron interactions are supposed to be strictly local. There, the Wannier functions [Equation (7)] would be very different for the occupied and the unoccupied states, especially for a charge-transfer gap insulator where the occupied dn-1 state is inside the O—2p band while the unoccupied dn+1 state is very well above. In this situation, there would be no easy way to define properly the Hubbard interaction parameter U.
As an example of how much spatial extension a nominally maximally localized atomic Wannier function can have and how strongly it can be affected by covalence effects, let us consider Ni—3d orbitals in the charge-transfer gap insulator NiO. Figure 3 shows the NiO minority spin band structure calculated for ferromagnetically aligned Ni magnetic moments using LDA + U with JH = 1 eV and U = 2 eV [panel (a)], U = 12 eV [panel (b)], and U = 6 eV [panels (c)-(e)] applied to Ni—3d electrons. Fat bands indicate the strength of the Ni—3d character, while the dashed lines are the eigenvalues of an effective Hamiltonian in the basis of maximally localized Ni—t2g Wannier functions obtained using the Wannier90 package . For the very small U value of 2 eV and the very large U value of 12 eV, the Ni—t2g bands lie either above or below the oxygen—2p bands, respectively, which reduces their hybridization and makes it possible to obtain a reasonable effective Hamiltonian with eigenvalues closely matching the LDA + U bands. However, the resulting Wannier functions in real space (shown below the corresponding band-structure plots in Figure 3) look quite different in the two cases and, in particular, have different degrees of spatial extension. In the most realistic case of U = 6 eV, the Ni—t2g and oxygen—2p bands are
Figure 3. Maximally localized Ni—t2g Wannier functions in NiO. At the top shown are the NiO spin minority band structures, with fat bands indicating the Ni—3d orbital character, as well as the eigenvalues of the Ni—t2g maximally localized Wannier functions based effective Hamiltonians. The Fermi level is set to zero and marked with a dashed black line. One of the Ni—t2g Wannier functions in real space is shown below each corresponding band-structure plot. The band-structure results were obtained in LDA + U with (a) U = 2 eV, (b) U = 12 eV, and (c)-(e) U = 6 eV. In (c)-(e), the effective Hamiltonian’s eigenvalues were constrained during wannierization to match (c) the lowest, (d) the highest, or (e) no LDA + U bands in the occupied manifold.
energetically very close and so get strongly hybridized. As a result, wannierization may give even more ambiguous results in this case. Depending on whether certain constraints are applied [Figure 3(c) and Figure 3(d)] or not [Figure 3(e)], one can obtain either very extended and ill-shaped or very well-localized Wannier functions. In the latter case, however, the effective Hamiltonian’s eigenvalues are a poor match to the LDA + U bands.
The cuprates are another good example where we can define the upper Hubbard band very well because it is well separated from the oxygen—2p band, while the lower Hubbard band, involving also all kinds of multiplets spread out over about 8 eV, is for a large part inside the oxygen—2p band. Also, the lowest energy d8 state without the hybridization switched on is a triplet F state and not a singlet state, as one assumes in a Hubbard model. In fact, the singlet A1gd8 state hybridizes so strongly that it pushes out the Zhang-Rice singlet state from the top of the O band and this state has more oxygen—2p than Cu—3d character in it. This is very similar to the BaBiO3 “bound two-hole state” pushed up above the Fermi level that we discussed in Section 3. If the Zhang-Rice state is well enough pushed out of the top of the O band, then we could actually define a single site Wannier function for this but it would have more density on O than on Cu.
In this case of correlated oxides where transition metal—d orbitals are strongly localized yet also subject to hybridization with oxygen—p orbitals, a way to improve on the localization of Wannier functions would be to also include oxygen—2p orbitals into the Wannier basis. This, however, may significantly increase the size of the Hilbert space required in model calculations. On the other hand, the increase can be much less dramatic if only the most important molecular combinations of oxygen—p orbitals are considered, such as the oxygen—(x2-y2) or oxygen—a1g in our earlier examples of NaNiO2 and BaBiO3. Recently, this promising MO based approach has been successfully applied to calculate resonant X-ray spectral responses in rare-earth nickelates . Another useful application would be to study systems with strong electron-phonon coupling, especially of the kind that strongly affects hybridization between cation and oxygen orbitals. There is, for example, a strong electron-phonon coupling of this kind in BaBiO3 where the A1g—symmetric (so-called “breathing”) oxygen phonon mode is coupled to the hybridization strength between the Bi-6s atomic orbital and the oxygen—a1g MO. The role of this coupling in the bismuthates’ superconductivity has recently been a subject of intense theoretical research, both at the band theory    and also model Hamiltonians levels, the latter including conventional atomic orbital based models  as well as MO based ones .
5. Strongly Correlated Molecular Orbitals (Na2IrO3 as an Example)
Since molecular orbitals are localized objects, electrons or holes occupying them may be subject to strong on-site Coulomb repulsion, where a site can now be a molecular orbital one and technically comprise more than one atomic site. From the band theory point of view, it would be very desirable to have a mean-field method, perhaps similar to the conventional LDA + U, in order to describe local correlations in systems with dominant MO character of the valence bands. Interestingly, application of the conventional LDA + U method to such systems may have an unphysically detrimental effect on their MO nature.
As an illustration, let us consider sodium iridium oxide Na2IrO3, a famous candidate Kitaev system , with an unusual property that its Ir—5d orbitals hybridize in a way such as to form MOs on Ir hexagons  . It is worth emphasizing that, in contrast with our previous examples, here MOs are formed out of cation instead of oxygen orbitals, which makes this case rather special. This becomes possible thanks to the larger spatial extent of the 5d orbitals, compared to that of the 3d ones, which leads to their larger nearest-neighbor hybridization, and also thanks to the peculiar geometrical arrangement of the IrO6 octahedra. Figure 4(a) shows the Na2IrO3 DOS projected onto Ir—5d MOs, as obtained from a generalized gradient approximation  (GGA) calculation. Here, we would like to focus on the MO aspects of the Na2IrO3 electronic structure and therefore neglect spin-orbit coupling effects, which one would otherwise need to also take into account. In this calculation, we assume a zig-zag order of Ir magnetic moments, shown in Figure 4(b), since it is experimentally found to be the ground magnetic state for Na2IrO3 . As one can see from the GGA calculations, this magnetic order can be well described in terms of on-hexagon molecular orbitals. In particular, near the Fermi level we find a completely filled A1g + E2u MO and slightly filled A1g − E2u MO in one spin channel and an opposite situation in the other spin channel, with A1g ± E2u standing for linear combinations of MOs with the A1g and E2g symmetries. The connection between the zig-zag magnetic order and such peculiar MO occupations becomes
Figure 4. Ir molecular orbitals in Na2IrO3. (a) and (e) show Ir MO projected DOS of Na2IrO3 calculated using GGA and GGA + U, respectively. (b) The zig-zag order of Ir magnetic moments within a hexagonal Ir layer. (c) and (d) show, respectively, the isolated A1g + E2u and A1g − E2u Ir MOs.
clear if one inspects the shapes of these two MOs in real space. Indeed, as shown in Figure 4(c) and Figure 4(d), one of these MOs is located on the right side of an Ir hexagon while the other MO is located on its left side, and this exactly produces the zig-zag order if the two MOs are occupied by holes of opposite spins. Now, there is one more MO at the Fermi level, namely, E2u, which is slightly empty in both spin channels. If a method like GGA + U but designed to act in the basis of MOs were applied to Na2IrO3, it would push up in energy and completely empty the A1g − E2u and A1g + E2u MO states and would push down and complete fill the E2g MO states. This is not what happens when the conventional GGA + U method (U = 2.7 eV and JH = 0.7 eV on Ir-5d orbitals) is applied [Figure 4(e)]. Although it does open a charge gap of 0.25 eV, but in a way that splits and mixes the A1g − E2u, A1g + E2u, and E2u MO states, thus destroying the MO nature of the Na2IrO3 valence bands. Similar damage is also done to the lower-lying B1u and E1g MOs.
For the local MO correlations to be treated properly (at least, in a mean-field sense), two approaches seem to be in place. In the first approach, one can first construct MOs via the unitary basis transformation of Equation (2) and then use the conventional LDA + U expressions for total energy and potential but written in the basis of MOs. We see three problems associated with this approach. First, it is not always obvious in which way MOs should be constructed. Second, this approach is intrinsically rotationally non-invariant, i.e., would produce different results depending on the choice of the basis. Third, calculation of the Coulomb and exchange interaction matrix elements and , with Mi and Mj denoting MOs, may pose a serious challenge.
The second and, as we believe, more promising approach consists in extending the atomic orbital LDA + U method by additionally taking into account Coulomb interactions between different atomic sites. Also referred to as extended DFT + U + V, this approach was first discussed by Campo and Cococcioni in . It can be shown that if all important inter-site interaction terms are taken into account, this approach can capture the same MO physics as the more intuitive first approach discussed above. Whether the inter-site interaction terms that the authors of  chose to keep in their interaction Hamiltonian are sufficient in doing so is an open question, and further investigations into this matter would be of great value. By the way, recognizing the importance of inter-site interactions explains why GGA + U fails in the case of Na2IrO3. Indeed, although GGA suffers from the self-interaction problem, it at least treats both on-site and inter-site interactions on the same footing. This allows GGA to capture qualitatively all MO energy splittings in Na2IrO3 [Figure 4(a)]. On the other hand, in GGA + U the on-site interactions are unphysically exaggerated in comparison with the inter-site ones, which apparently must be similarly strong.
In summary, we presented our band theory perspective on some important aspects of molecular orbitals in complex oxides. In particular, we first discussed technical issues related with the implementation of MO-based projection of electronic states, using NaNiO2 as an example. Then, this discussion was extended by considering non-orthogonal MOs, as observed in BaBiO3, and the way of performing their orthonormalization. This led us to the concept of atomic and molecular orbital Wannier functions, where we discussed the problem of the Wannier functions being typically rather extended in real space and how this property may restrict their usage as basis functions in model Hamiltonians, using NiO as an example. A suggestion was then made that MO based model Hamiltonians might be a good approach to describing correlated transition metal systems with strong cation-anion orbital hybridization and also systems with strong electron-phonon coupling. Our final point concerned the problem of a proper band theory treatment of strong correlations between electrons or holes occupying molecular rather than atomic orbitals. In this context, we discussed the failure of the conventional atomic LDA + U method to capture the MO nature of the electronic states in Na2IrO3 and outlined possible pathways towards a more adequate band theory based description. As the interest in better understanding of the role of molecular orbitals in complex oxides is growing, we expect that our current review will provide a useful reference for future studies.
This work was supported by NSERC, CIfAR, and the Max Planck—UBC Stewart Blusson Quantum Matter Institute.