A Band Theory Perspective on Molecular Orbitals in Complex Oxides

Kateryna Foyevtsova^1,2, George A. Sawatzky^1,2

¹Department of Physics & Astronomy, University of British Columbia, Vancouver, BC, Canada.
²Stewart Blusson Quantum Matter Institute, University of British Columbia, Vancouver, BC, Canada.

DOI: 10.4236/jmp.2019.108062                

Abstract

In view of the growing interest in molecular orbitals (MOs) encountered in certain complex oxides, we review some of their properties from the band theory perspective and provide detailed examples based on real materials. Our discussion includes some technical aspects of identifying MOs in electronic structure calculations and considers cases when MOs can be both orthogonal and non-orthogonal. We also describe orthonormalization of MOs, a procedure converting them into Wannier functions, and discuss the problem of Wannier functions possibly being rather spatially extended and how using MO, rather than atomic orbital, based effective Hamiltonians might be a better choice in describing certain strongly correlated systems as well as systems with strong electron-phonon coupling. Furthermore, we address the problem of strongly correlated MOs and how it can be treated in band theory calculations.

Keywords

Electronic Structure of Solids, Molecular Orbitals, Complex Oxides, Negative Charge-Transfer Insulators

1. Introduction

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 [1]. 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 [2] [3] and superconducting bismuthates [4] [5] [6]. 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 [7]. Among this important class of materials are such vigorously discussed but still controversial systems as superconducting cuprates [8] [9] and bismuthates [10] [11] [12], rare-earth nickelates [13], and transition-metal oxide based Li-ion battery cathode materials [14] [15]. 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 (NaNiO₂ as an Example)

Let us consider the prototypical negative charge-transfer insulator NaNiO₂ as our first example. Its monoclinic C2/m crystal structure consists of triangular layers of edge-sharing cooperatively elongated NiO₆ octahedra intercalated with Na⁺ ions, with one formula unit per primitive unit cell [16]. 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 Ni²⁺L, where L is an oxygen hole [17]. This three-hole state has a net spin of ½ and would formally correspond to Ni³⁺ ( $e_g^3$ ) with spin 1/2, which is quite contrary to what one would expect from Hund’s rules if it really were Ni³⁺. Due to the strong hybridization between the Ni—e_g orbitals and the oxygen MOs of respective symmetry, the oxygen holes select to occupy the O—(x²-y²) MOs, one hole per oxygen octahedron, and form a Zhang-Rice-like spin singlet state with the hole in the Ni—(x²-y²) orbital. This selective occupation of the (x²-y²) MOs by the holes, which results in the mentioned cooperative elongation of the NiO₆ 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 NaNiO₂ 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 $\frac{1}{4}{J}_{pd}$ rather than $\frac{3}{4}{J}_{pd}$ , where J_pd 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 NaNiO₂ 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 NaNiO₂. (a) Projected densities of states (PDOS) per formula unit for the Ni—3d atomic orbitals and the oxygen MOs. (b) An isolated oxygen—(x²-y²) MO ${\psi }_R^M\left(r\right)$ . 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—(x²-y²) character highlighted in red. (d) An oxygen—(x²-y²) Bloch function ${\psi }_R^M\left(r\right)$ at the $\Gamma$ 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 $|p_x\rangle$ , $|p_y\rangle$ , and $|p_z\rangle$ 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) [18] [19] calculation (on-site Coulomb repulsion U = 6 eV, Hund’s exchange interaction J_H = 1 eV for Ni—3d electrons [20] ) performed, like the rest of the calculations presented in this paper, using the linearized augmented plane waves method implemented in the Wien2k package [21]. 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—(x²-y²) atomic orbital and the oxygen—(x²-y²) 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—(x²-y²) molecular orbital ${\psi }_R^M\left(r\right)$ 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 ${\psi }_k^M\left(r\right)$ constructed from ${\psi }_R^M\left(r\right)$ :

${\psi }_k^M\left(r\right)={\displaystyle {\sum }_R{\psi }_R^M\left(r\right){\text{e}}^{-ikR}}$ , (1)

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—(x²-y²) MO Bloch function at the $\Gamma$ point [k = (0, 0, 0)] in Figure 1(d).

Generally speaking, MO Bloch functions ${\psi }_k^{M_i}\left(r\right)$ can be obtained from the atomic orbital Bloch functions ${\psi }_k^{a_j}\left(r\right)$ via a unitary transformation U(k):

${\psi }_k^{M_i}\left(r\right)={\displaystyle {\sum }_j{U}_{ij}\left(k\right){\psi }_k^{a_j}\left(r\right)}$ . (2)

Allowing for the k-dependence of the transformation matrix elements U_ij(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—(x²-y²) 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 NiO₆ octahedron are in fact the $|p_x\rangle$ , $|p_y\rangle$ , and $|p_z\rangle$ 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 NaNiO₂.

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 NaNiO₂, where indeed the majority spin hole state has a strong oxygen—(x²-y²) 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—(x²-y²) MO and the Ni—(x²-y²) atomic orbital in NaNiO₂], there are a number of notable exceptions. For example, in the iron disulfide FeS₂ (pyrite structure) and the recently discussed iron dioxide FeO₂ [22], 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 KO₂, 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 $|u_l^{\alpha ,\sigma }\left(E_{1l}\right)Y_m^l\rangle$ , which are the solutions of the Schrödinger equation within the muffin-tin sphere of atom α at the linearization energy $E_{1l}$ and where $|Y_m^l\rangle$ 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 $|Y_m^l\rangle$ inside muffin-tin spheres. In practice, the two approaches give very similar results for the projected DOS, but having also the radial part $u_l^{\alpha ,\sigma }\left(E_{1l},r\right)$ 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 (BaBiO₃ as an Example)

For certain lattice structures—including those of the cuprates, with their characteristic two-dimensional square CuO₂ 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 BaBiO₃. In a perovskite lattice structure of the ABO₃ type [Figure 2(a)], neighboring BO₆ 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 BaBiO₃, the MO of most importance is oxygen—a_1g, 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—a_1g 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, $\langle {\psi }_k^M|{\psi }_k^M\rangle$ , 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—a_1g MO in BaBiO₃ at the $\Gamma$ 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 [1], the k-vector dependent normalization factor ${\beta }_k$ ,

${\left|{\beta }_k\right|}^2\langle {\psi }_k^M|{\psi }_k^M\rangle =1$ , (3)

Figure 2. Oxygen molecular orbitals in BaBiO₃. (a) An idealized (e.g., without octahedra’s rotations or bond disproportionation) cubic perovskite crystal structure of BaBiO₃, with one formula unit per unit cell, featuring also an isolated oxygen—a_1g MO ${\psi }_R^M\left(r\right)$ . (b) Projection of the BaBiO₃ electronic states onto oxygen MOs; the Fermi level is marked with a dashed line. (c) and (d) show oxygen—a_1g MO Bloch functions at the $\Gamma$ [k = (0, 0, 0)] and R [k = (π, π, π)] points, respectively. Note that at the $\Gamma$ point the individual ${\psi }_R^M\left(r\right)$ contributions completely cancel out due to their non-orthogonality, which results in vanishing oxygen—a_1g MO weight at this k-vector in (b).

may be ill-defined at certain k-points. For the oxygen—a_1g MO in BaBiO₃, for example, it diverges at the $\Gamma$ point:

${\beta }_k={\left[1-\frac{1}{3}\left(\cos k_x+\cos k_y+\cos k_z\right)\right]}^{-1/2}$ , (4)

as obtained using the oxygen—a_1g MO given by

$|{\psi }_k^{a_{1g}}\rangle =\frac{1}{\sqrt{6}}\left[|p_k^x\rangle \left(1-{\text{e}}^{-ia{k}_x}\right)+|p_k^y\rangle \left(1-{\text{e}}^{-ia{k}_y}\right)+|p_k^z\rangle \left(1-{\text{e}}^{-ia{k}_z}\right)\right]$ , (5)

where $|p_k^y\rangle$ , 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—a_1g MO and the Bi—6s atomic orbital $|s_k\rangle$ ,

$|{\psi }_k^{s-a_{1g}}\rangle =\frac{1}{\sqrt{{\alpha }^2+{\beta }^2}}\left(\alpha |s_k\rangle \pm \beta |{\psi }_k^{a_{1g}}\rangle \right)$ , (6)

which was used by us and our co-authors in [5] for analyzing the BaBiO₃ 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 ${\psi }_R^M\left(r\right)$ of any kind via multiplication of its corresponding Bloch function ${\psi }_k^M\left(r\right)$ by ${\beta }_k$ converts this molecular orbital into a Wannier function [23]. 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 $|{\Psi }_{k,\nu }\rangle$ :

$|w_k\rangle ={\displaystyle {\sum }_{n<\nu <m}\langle {\Psi }_{k,\nu }|{\psi }_k\rangle |{\Psi }_{k,\nu }\rangle }$ . (7)

Here, $|{\psi }_k\rangle$ can be either an atomic or a molecular orbital Bloch function and $\nu$ is a band index running from band n to band m. On top of this construction, one may also apply the maximal localization [23] and the disentanglement [24] procedures. This is finalized by orthonormalization [Equation (3)] to determine ${\beta }_k$ . The resulting Wannier functions $|w_k\rangle$ 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 ${\beta }_k$ 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 $|w_k\rangle$ [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 $|w_k\rangle$ [Equation (7)] would be very different for the occupied and the unoccupied states, especially for a charge-transfer gap insulator where the occupied d^n-1 state is inside the O—2p band while the unoccupied dⁿ⁺¹ 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 $|w_k\rangle$ 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 J_H = 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—t_2g Wannier functions obtained using the Wannier90 package [25]. For the very small U value of 2 eV and the very large U value of 12 eV, the Ni—t_2g 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—t_2g and oxygen—2p bands are

Figure 3. Maximally localized Ni—t_2g 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—t_2g 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—t_2g 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 d⁸ 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 A_1gd⁸ 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 BaBiO₃ “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—(x²-y²) or oxygen—a_1g in our earlier examples of NaNiO₂ and BaBiO₃. Recently, this promising MO based approach has been successfully applied to calculate resonant X-ray spectral responses in rare-earth nickelates [3]. 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 BaBiO₃ where the A_1g—symmetric (so-called “breathing”) oxygen phonon mode is coupled to the hybridization strength between the Bi-6s atomic orbital and the oxygen—a_1g MO. The role of this coupling in the bismuthates’ superconductivity has recently been a subject of intense theoretical research, both at the band theory [4] [5] [6] and also model Hamiltonians levels, the latter including conventional atomic orbital based models [26] as well as MO based ones [27].

5. Strongly Correlated Molecular Orbitals (Na₂IrO₃ 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 Na₂IrO₃, a famous candidate Kitaev system [28], with an unusual property that its Ir—5d orbitals hybridize in a way such as to form MOs on Ir hexagons [29] [30]. 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 IrO₆ octahedra. Figure 4(a) shows the Na₂IrO₃ DOS projected onto Ir—5d MOs, as obtained from a generalized gradient approximation [31] (GGA) calculation. Here, we would like to focus on the MO aspects of the Na₂IrO₃ 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 Na₂IrO₃ [32]. 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 A_1g + E_2u MO and slightly filled A_1g − E_2u MO in one spin channel and an opposite situation in the other spin channel, with A_1g ± E_2u standing for linear combinations of MOs with the A_1g and E_2g symmetries. The connection between the zig-zag magnetic order and such peculiar MO occupations becomes

Figure 4. Ir molecular orbitals in Na₂IrO₃. (a) and (e) show Ir MO projected DOS of Na₂IrO₃ 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 A_1g + E_2u and A_1g − E_2u 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, E_2u, 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 Na₂IrO₃, it would push up in energy and completely empty the A_1g − E_2u and A_1g + E_2u MO states and would push down and complete fill the E_2g MO states. This is not what happens when the conventional GGA + U method (U = 2.7 eV and J_H = 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 A_1g − E_2u, A_1g + E_2u, and E_2u MO states, thus destroying the MO nature of the Na₂IrO₃ valence bands. Similar damage is also done to the lower-lying B_1u and E_1g 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 $U_{M_i{M}_j}$ and $J_{M_i{M}_j}$ , with M_i and M_j 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 [33]. 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 [33] 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 Na₂IrO₃. 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 Na₂IrO₃ [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.

6. Conclusion

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 NaNiO₂ as an example. Then, this discussion was extended by considering non-orthogonal MOs, as observed in BaBiO₃, 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 Na₂IrO₃ 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.

Acknowledgements

This work was supported by NSERC, CIfAR, and the Max Planck—UBC Stewart Blusson Quantum Matter Institute.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Zhang, F.C. and Rice, T.M. (1988) Physical Review B: Condensed Matter, 37, 3759-3761. https://doi.org/10.1103/PhysRevB.37.3759
[2]	Johnston, S., Mukherjee, A., Elfimov, I., Berciu, M. and Sawatzky, G.A. (2014) Physical Review Letters, 112, Article ID: 106404. https://doi.org/10.1103/PhysRevLett.112.106404
[3]	Green, R.J., Haverkort, M.W. and Sawatzky, G.A. (2016) Physical Review B: Condensed Matter, 94, Article ID: 195127. https://doi.org/10.1103/PhysRevB.94.195127
[4]	Foyevtsova, K., Khazraie, A., Elfimov, I. and Sawatzky, G.A. (2015) Physical Review B: Condensed Matter, 91, Article ID: 121114. https://doi.org/10.1103/PhysRevB.91.121114
[5]	Khazraie, A., Foyevtsova, K., Elfimov, I. and Sawatzky, G.A. (2018) Physical Review B: Condensed Matter, 97, Article ID: 075103. https://doi.org/10.1103/PhysRevB.97.075103
[6]	Khazraie, A., Foyevtsova, K., Elfimov, I. and Sawatzky, G.A. (2018) Physical Review B: Condensed Matter, 98, Article ID: 205104. https://doi.org/10.1103/PhysRevB.98.205104
[7]	Zaanen, J., Sawatzky, G.A. and Allen, J.W. (1985) Physical Review Letters, 55, 418-421. https://doi.org/10.1103/PhysRevLett.55.418
[8]	Bednorz, J.G. and Müller, K.A. (1986) Zeitschrift für Physik B, 64, 189-193. https://doi.org/10.1007/978-94-011-1622-0_32
[9]	Leggett, A.J. (2006) Nature Physics, 2, 134-136. https://doi.org/10.1038/nphys254
[10]	Cava, R.J., et al. (1988) Nature, 332, 814-816. https://doi.org/10.1038/332814a0
[11]	Sleight, A.W., Gillson, J.L. and Bierstedt, P.E. (1975) Solid State Communications, 17, 27-28. https://doi.org/10.1016/0038-1098(75)90327-0
[12]	Sleight, A.W. (2015) Physica C: Superconductivity and Its Applications, 514, 152-165. https://doi.org/10.1016/j.physc.2015.02.012
[13]	Medarde, M.L. (1997) Journal of Physics: Condensed Matter, 9, 1679-1707. https://doi.org/10.1088/0953-8984/9/8/003
[14]	Liu, W., et al. (2015) Angewandte Chemie International Edition in English, 54, 4440-4457. https://doi.org/10.1002/anie.201409262
[15]	Okubo, M. and Yamada, A. (2017) ACS Applied Materials & Interfaces, 9, 36463-36472. https://doi.org/10.1021/acsami.7b09835
[16]	Dyer, L.D., Borie, B.S. and Pedro Smith, G. (1954) Journal of the American Chemical Society, 76, 1499-1503. https://doi.org/10.1021/ja01635a012
[17]	Kuiper, P., Kruizinga, G., Ghijsen, J., Sawatzky, G.A. and Verweij, H. (1989) Physical Review Letters, 62, 221-224. https://doi.org/10.1103/PhysRevLett.62.221
[18]	Perdew, J.P. and Wang, Y. (1992) Physical Review B: Condensed Matter, 45, 13244-13249. https://doi.org/10.1103/PhysRevB.45.13244
[19]	Liechtenstein, A.I., Anisimov, V.I. and Zaanen, J. (1995) Physical Review B: Condensed Matter, 52, R5467-R5470. https://doi.org/10.1103/PhysRevB.52.R5467
[20]	Anisimov, V.I., Zaanen, J. and Andersen, O.K. (1991) Physical Review B: Condensed Matter, 44, 943-954. https://doi.org/10.1103/PhysRevB.44.943
[21]	Schwarz, K., Blaha, P. and Madsen, G.K.H. (2002) Computer Physics Communications, 147, 71-76. https://doi.org/10.1016/S0010-4655(02)00206-0
[22]	Streltsov, S.S., Shorikov, A.O., Skornyakov, S.L., Poteryaev, A.I. and Khomskii, D.I. (2017) Scientific Reports, 7, Article No. 13005. https://doi.org/10.1038/s41598-017-13312-4
[23]	Marzari, N. and Vanderbilt, D. (1997) Physical Review B, 56, 12847-12865. https://doi.org/10.1103/PhysRevB.56.12847
[24]	Souza, I., Marzari, N. and Vanderbilt, D. (2001) Physical Review B: Condensed Matter, 65, 9028. https://doi.org/10.1103/PhysRevB.65.035109
[25]	Mostofi, A.A., et al. (2014) Computer Physics Communications, 185, 2309-2310. https://doi.org/10.1016/j.cpc.2014.05.003
[26]	Li, S. and Johnston, S. (2019) Quantum Monte Carlo Study of Lattice Polarons in the Two-Dimensional Multi-Orbital Su-Schrieffer-Heeger Model.
[27]	Yam, Y.-C., Berciu, M. and Sawatzky, G.A. to be published.
[28]	Chaloupka, J., Jackeli, G. and Khaliullin, G. (2010) Physical Review Letters, 105, Article ID: 027204. https://doi.org/10.1103/PhysRevLett.105.027204
[29]	Mazin, I.I., Jeschke, H.O., Foyevtsova, K., Valentí, R. and Khomskii, D.I. (2012) Physical Review Letters, 109, Article ID: 197201. https://doi.org/10.1103/PhysRevLett.109.197201
[30]	Foyevtsova, K., Jeschke, H.O., Mazin, I.I., Khomskii, D.I. and Valentí, R. (2013) Physical Review B: Condensed Matter, 88, Article ID: 035107. https://doi.org/10.1103/PhysRevB.88.035107
[31]	Perdew, J.P., Burke, K. and Ernzerhof, M. (1996) Physical Review Letters, 77, 3865-3868. https://doi.org/10.1103/PhysRevLett.77.3865
[32]	Choi, S.K., et al. (2012) Physical Review Letters, 108, Article ID: 127204. https://doi.org/10.1103/PhysRevLett.108.127204
[33]	Campo, V.L. and Cococcioni, M. (2010) Journal of Physics: Condensed Matter, 22, Article ID: 055602. https://doi.org/10.1088/0953-8984/22/5/055602