Cooperative effects of strain and electron correlation in epitaxial VO$2$ and NbO$2$

Enormous efforts have been made to understand vanadium dioxide (VO$2$⁠), which exhibits a metal-to-insulator transition (MIT) near room temperature.¹ The MIT in VO$2$ is accompanied by a structural transition from the rutile R (metallic) to the monoclinic M$1$ (insulating) crystal structure, and it is understood that both Mott’s (electron correlation)^2,3 and Peierls’s (structural distortion)⁴ physics play essential roles in this transition. A large amount of theoretical research has been performed to investigate the electronic structures of VO$2$ using the local density approximation (LDA), LDA + U, LDA + DMFT, quantum Monte Carlo, etc. However, the role of electron correlations in driving the MIT is still under debate, due to the fact that the Mott insulating phase cannot be disentangled from the Peierls instability.^5–11 

Recent progress in strain-engineering VO$2$ grown by molecular beam epitaxy on TiO$2$ substrates has opened the possibility of modulating electron correlation with strain before the Peierls instability sets in Refs. 12–16. Since TiO$2$ remains in the rutile crystal structure at all temperatures and has a $c$-axis lattice constant approximately $3.8%$ longer than bulk VO$2$⁠, the strain can be engineered by choosing the thin film growth direction. If the growth direction is along the rutile $c$-axis, denoted as VO$2$ (001), the strain is small and the corresponding VO$2$ thin film behaves like bulk VO$2$⁠.¹⁷ On the other hand, if the growth direction is perpendicular to the rutile $c$-axis, e.g., (100) and (110), a large strain is induced due to significant elongation of the $c$-axis lattice constant. Recent evidence of a strain-induced orbital selective Mott state was reported from hard x-ray photoelectron spectroscopy (HAXPES) and x-ray absorption spectroscopy (XAS) experiments in epitaxial (100) and (110) films by Mukherjee et al.¹⁴ and Quackenbush et al.¹⁵ The reports revealed how electron correlation can be strain-enhanced in VO$2$ thin films, but a full theoretical description of the interplay between strain and electron correlation remained incomplete.

In this paper, we employ the density functional theory implemented with the slave spin method (DFT+SS) to investigate the electronic structure of VO$2$ thin films under strain. We focus on the rutile phase and study the phase diagram as a function of doping, strain, and correlation strength. We find that while the band structure is modified by the elongation of the $c$-axis driving more electrons into the $d∥$ band, the electron correlation has the opposite effect and tends to distribute electrons more evenly among the three $t2g$ bands. The interplay between these two effects results in orbital-dependent redistributions of the electron occupations and spectral weights, which pushes the strained VO$2$ toward an orbital selective Mott transition (OSMT). Our theoretical models are in excellent agreement with V L-edge resonant x-ray photoemission spectroscopy of both (001) and (100) epitaxial 10 nm VO$2$/TiO$2$⁠, which highlights the evolution of the lower Hubbard band as a function of strain. We have further employed the same theoretical approach to study rutile NbO$2$ under similar strain conditions and concluded that no OSMT can be found for comparable strain regimes. Our conclusion is supported by the high-temperature hard x-ray photoelectron spectroscopy (HAXPES) measurement on relaxed NbO$2$ (110) thin film with no trace of the lower Hubbard band. Our results demonstrate that electron correlation has significant effects on the physical properties of VO$2$ even in the rutile phase and can be modified dramatically by strain engineering.

The Mott physics can be generally described by a Hamiltonian containing two terms: $H=Ht+HU$⁠, where $Ht$ is the kinetic energy, typically obtainable from band structure calculations, and $HU$ is the multi-orbital Hubbard interaction that provides a Coulomb penalty for two electrons on the same atom. In the single band case, if the system is at “half-filling,” which refers to the electron density being 1 electron per site (⁠$n=1$⁠), the transition from metallic to the Mott insulating state can be tuned by varying the ratio of $U/W$⁠, where $U$ is the Hubbard on-site Coulomb energy and $W$ is the bandwidth.

The physics of “Mottness” is even richer in a multi-orbital system. The half-filling condition for the Mott insulating state becomes $n=No$⁠, where $No$ is the number of orbitals, but there could exist a new state beyond the standard Mott insulator. Let us consider a generic two-orbital system with $dxz$ and $dyz$⁠, whose typical Fermi surface is plotted in Fig. 1(a). At half-filling (⁠$n=2$⁠), the Mott insulating state can occur with a large enough on-site Coulomb interaction, and the Fermi surface will be completely gapped out, as shown in Fig. 1(b). As the electron density moves away from half-filling, the correlation effect becomes weaker, and it is generally expected that the Mott physics becomes much less important. However, strain engineering facilitates the OSMT mechanism as explained below.

Let us assume that $nxz,yz$ and $Wxz,yz$ are the electron density and bandwidth of orbitals $dxz$ and $dyz$⁠, respectively. The total electron density is just $n=nxz+nyz$⁠, and the correlation effect on each orbital can be characterized by the ratio of $U/Wxz,yz$⁠. The key part of this mechanism is that the orbital wavefunctions are highly anisotropic in real space, thus the bandwidths can be modulated quite differently in different orbitals under the same tensile strain. In this example of the two-orbital system, the $dxz$ (⁠$dyz$⁠) orbital has a much larger wavefunction overlap along the $x$ (⁠$y$⁠) direction. If a tensile strain is applied along the $x$ direction, the lattice constant along $x$ direction will be elongated, which reduces $Wxz$ more significantly than $Wyz$⁠. This orbital-dependent bandwidth reduction due to tensile strain leads to two important consequences. First, it is entirely possible to put more electrons in the $dxz$ orbital. Second, the correlation effect is increased only in the $dxz$ orbital, expressed by the ratio of $U/Wxz$⁠. As a result, even if the system is not strictly half-filled originally, the application of the strain along the $x$ direction could drive the system to a state where only the $dxz$ orbital is Mott insulating (⁠$nxz∼1$ with larger $U/Wxz$⁠), while the $dyz$ orbital remains weakly-correlated (⁠$nyz$ away from $1$ with smaller $U/Wyz$⁠). In this new state of matter, namely, the orbital selective Mott state, the Fermi surfaces will be disconnected, exhibiting an interesting “Fermi arc” profile shown in Fig. 1(c). The mechanism described above is quite general for a system with multiple orbitals near the Fermi surface together with strong Hubbard interactions.^18–23 In Sec. III, we present our theoretical study and supporting experimental data for the case of VO$2$ thin films. We also perform the same calculations for rutile NbO$2$ under similar strain conditions. Our results show that the OSMT occurs in VO$2$ thin films, but not in NbO$2$⁠, which can be attributed to the fact that NbO$2$ has much larger bandwidths in general compared to VO$2$⁠.

We solve the Hamiltonian expressed as

Here, $HTB$ is the tight-binding model composed of the $d$-orbtials of the V atom as well as the $p$ orbitals of the O atom, and the hopping parameters are determined by fitting the band structure calculated by the density functional theory using Wannier90.²⁴ The DFT calculations are performed with full potential linear augmented plane waves plus local orbitals (FP-LAPW+lo) and the Perdew-Burke-Ernzerhof generalized gradient approximation (PBE-GGA) provided in the WIEN2k code.²⁵ The electron correlations are treated by the multi-orbital Hubbard term $HU$ on the $d$ orbitals of the V atom defined as

where $ciασ†$ creates an electron in a V $d$ orbital $α$ with spin $σ$ at site $i$⁠. $J$ is Hund’s coupling, $U′$ the inter-orbital coupling, $σ¯=−σ$⁠, $n^iασ=ciασ†ciασ$⁠, and $Hpair$ is the pair hopping term defined as

Lastly,

is the on-site energy introduced to make $HU=0$ at the half-filling condition.²⁶ It is easy to show that for a five orbital model, $μdU=(9U−20J)/4$⁠, if we adopt the form of the Hubbard interactions defined in Eq. (2). We will use the values of $U=6$ eV, $J=1$ eV, and $U′=U−2J$ for VO$2$ and NbO$2$ throughout this paper.

The slave spin method^{14,23,26–38} will be employed to treat $H=HTB+HU$⁠. It is worth mentioning that the slave spin formalism^23,28 has been shown to capture the featureless Mott transition in good agreement with DMFT results, and it correctly reproduces the Gutzwiller limit for the quasiparticle weight $Z→2x/(1+x)$⁠, where $x$ is the doping away from the half-filling in the large $U$ limit.^14,28,31 We adopt the $U(1)$ version of the slave spin method described in Ref. 31 and focus only on the rutile phase of VO$2$ in the normal state without any spontaneous symmetry-breaking order in this study. A brief review on the formalism of $U(1)$ slave spin method can be found in  Appendix A.

The strained case of VO$2$ thin films grown on a TiO$2$ substrate with the growth direction along the rutile $a$-axis, denoted as VO$2$ (100), is schematically illustrated in Fig. 2. In this case, the lattice constants of the $b$- and $c$-axes of the VO$2$ match those of the TiO$2$ substrate, while the $a$-axis lattice constant shrinks. Table I summarizes the lattice constants used in the present study. We now compare the unstrained (bulk) VO$2$ and VO$2$ (100) to resolve the interplay between band structures and electron correlations.

We first present our results without the interaction term $HU$⁠. This will serve as a reference to highlight the effects of strain on the band structures. We follow the local coordinate system introduced by Eyert,⁵ in which the $t2g$ orbitals are $dx2−y2$ (⁠$d∥$⁠), $dxz$⁠, and $dyz$ (⁠$dπ$ bands) and the $eg$ orbitals are $d3z2−r2$ and $dxy$⁠. Figure 3(a) summarizes our results for the density of states (DOS) from $HTB$ in unstrained VO$2$⁠. This result reproduces the previous DFT calculations accurately.⁵ The electron density in the vanadium $d$ and oxygen $p$ orbitals has been found to deviate significantly from the ionic picture, one that predicts full occupation of the oxygen $p$ orbitals, and 1 electron in the vanadium $d$ manifold. While DFT calculations obtain a total electron filling in V $d$ and O $p$ orbitals per unit cell (⁠$ntot$⁠) to be thirteen as expected, the electron filling in the vanadium $d$ orbitals (⁠$nd$⁠) is significantly larger than 1 due to strong hybridization between V $d$ and O $p$ orbitals. Previous calculations have obtained $nd$ ranging from $1.0$ to $3.15$⁠.^5,6,8,9 Because in DFT-SS we employ a tight-binding model, and the total electron filling can be varied freely by changing the chemical potential $μ$⁠, we can access general trends in the range of $13<ntot<14$⁠, which covers the range of $nd$ of the experimental interest.

Figures 3(c) and 3(d) plot the electron filling in each $d$ orbital as a function of $ntot$ for the unstrained VO$2$ and VO$2$ (100). Clearly, in both cases, the electron filling in the $eg$ orbitals remains almost the same as $ntot$ increases, while those in $t2g$ orbitals increases significantly as $ntot$ increases. This indicates that only the $t2g$ orbitals are important near the Fermi energy and consequently, the correlation effects will be primarily on $t2g$ orbitals.

Now, we discuss the effects of the (100) strain on the band structure. Comparing Figs. 3(a) and 3(b), it can be seen that under this strain, the bandwidth of the $dx2−y2$ (⁠$d∥$⁠) near the Fermi energy is significantly reduced, while those of the $dxz$ and $dyz$ (⁠$dπ$ bands) remain roughly the same. This orbital-dependent reduction of the bandwidth is one of the crucial ingredients for the occurrence of the orbital selective Mott state driven by strain.

Before we turn on the Hubbard interactions $HU$⁠, it is helpful to discuss the general effects of $HU$⁠. First, it costs energy to put two electrons on the same site, so electron hopping is suppressed. This results in the reduction of the quasiparticle weight $Z$⁠. Second, $HU$ could produce orbital-dependent effective on-site potentials that lead to the redistribution of orbital occupation numbers. To see why this happens, we analyze $HU$ given in Eq. (2). If $J$ is nonzero, $U′=U−2J<U$⁠. This means that if two electrons are at the same site, the interaction energy will be lower for them to stay in different orbitals than in the same orbital. As a result, $HU$ favors an equal orbital filling and this tendency is reflected in the orbital-dependent effective on-site potentials generated by $HU$⁠. Such an orbital-dependent on-site potential can be ignored only if the Hund’s coupling $J$ is negligible, which does not seem to be the case for VO$2$⁠.^8,9,11,39

The slave-spin method can capture both effects from $HU$⁠, as demonstrated in previous works.^{14,31,32,37,38} The orbital occupation numbers and the quasiparticle weights as a function of the total electron filling $ntot$ with $U=6$ eV and $J=1$ eV for the unstrained VO$2$ are plotted in Figs. 4(a) and 4(c). It can be seen that the Hubbard interactions make the electron occupation of the $t2g$ bands more symmetrical. Moreover, the quasiparticle weights of $t2g$ bands are much smaller than those of the $eg$ bands. These observations are consistent with the discussions given above. It is also interesting to note that the correlation effects are more and more pronounced as $ntot$ approaches 14. As mentioned above, since the $eg$ bands are away from the Fermi energy, they are not affected by $HU$⁠, and most extra electrons will go into the $t2g$ bands. Consequently, the electron filling in the $t2g$ bands increases with $ntot$⁠. If the electron filling in the $t2g$ bands (⁠$nt2g$⁠) approaches 3, the “half filling” condition for a three-orbital (⁠$t2g$⁠) model, the correlation effects in the $t2g$ bands will be enhanced significantly. From Fig. 4(a), it can be seen that $nt2g$ increases from 1.88 to 2.38 as $ntot$ increases from 13 to 14, which explains the rapid reduction of the quasiparticle weights in the $t2g$ bands.

Here, we discuss the strained case of VO$2$ (100) with Hubbard interactions. As mentioned above, the strain results in an orbital-dependent reduction of the bandwidth, and turning on the $HU$ interaction leads to even more non-trivial effects, as shown in Figs. 4(b) and 4(d). An intriguing difference in the behaviors of the quasiparticle weights can be seen from the comparison between strained and unstrained cases. As $ntot$ is fixed, the quasiparticle weights of the $dπ$ bands remain roughly the same, but the $d∥$ band decreases significantly in the strained system. This suggests that only the $d∥$ band is driven to be more correlated under strain. Moreover, we find that the electron filling in $d∥$ band increases, while that in $dπ$ bands decreases in VO$2$ (100). These two behaviors strongly indicate that strained VO$2$ is approaching an orbital selective Mott transition (OSMT), consistent with our previous reports^14,15

Figure 5 plots the differences in electron filling and quasiparticle weights between VO$2$ (100) and unstrained VO$2$⁠. In the non-interacting limit, shown in Fig. 5(a), we find that $Δn$ in the $d∥$ (red dots) changes sign as $ntot$ increases from 13 to 14. With the interactions, it can be seen in Figs. 5(b) and 5(c) that $Δn$ (⁠$ΔZ$⁠) in the $d∥$ increases (decreases) monotonically as $ntot$ increases from 13 to 14. This observation proves that electron correlation plays a crucial role here, and the mechanism for the OSMT relies on the cooperative interplay between strain-modulated band structures and electron correlation. In VO$2$ (100), because the $c$-axis is elongated due to the strain, and the $d∥$ band has a large wavefunction overlap along $c$-axis, the bandwidth of the $d∥$ band is significantly reduced compared to the $dπ$ bands. Since correlation effects are usually characterized by the ratio of the Hubbard interaction to the bandwidth (⁠$U/W$⁠), the $d∥$ band effectively becomes more correlated if Hubbard interactions exist. These two strain effects cooperatively induce the intriguing OSMT in VO$2$ (100).

Another important feature of the Mott physics is the emergence of the incoherent spectra, known as the upper and lower Hubbard bands (UHB and LHB), in the single particle spectral function.^40,41 We adopt the formalism of the single particle Green’s function in the slave-spin method derived in Ref. 30, and the spectral function at the mean-field level can be calculated directly by

where $Aα,σ(k→,ω)$ is the spectral function with orbital $α$ and physical spin $σ$⁠.³⁰ In the slave spin formalism, the electron Green’s function can be approximated as the convolution of the spinon and slave spin propagators, and contributions from the coherent and incoherent parts can be separately calculated.³⁰ The derivation of $Aα,σ(k→,ω)$ in the slave spin method can be found in  Appendix B. With the spectral function, we can simulate spectral functions for occupied states and unoccupied states by

Here, $nF(ω)$ is the Fermi Dirac function at room temperature $1β=25$ meV. The photoemission spectra is typically proportional to $Aoccu(ω)$⁠, the spectral weights of occupied states. Figure 6 presents the simulated spectra for cases of unstrained and VO$2$ (100) with several different total electron filling numbers. At smaller $ntot$⁠, the incoherent spectrum (green line) due to the LHB around 1 eV in energy can be clearly seen and becomes more pronounced in the strained sample. This is more evidence of strong correlation effects in the strained sample, consistent with the OSMT scenario. At larger $ntot$⁠, the incoherent spectrum is pushed to more negative energy, and both bulk and strained samples have a pronounced LHB. This is not surprising, since as $ntot$ approaches 14, the electron filling in the $t2g$ orbitals approach 3, which is “half-filling” for a three-orbital model. In this case, the system is already close to the Mott limit, even for the unstrained material. As a result, the physical properties do not change dramatically under the strain.

We previously presented evidence of a strain-induced OSMT for epitaxial VO$2$/TiO$2$ films in the high temperature, rutile phase using a combination of HAXPES and V L$3$-edge XAS.¹⁴ The preferential orbital filling was evidenced by the XAS dichroism, and the loss of quasiparticle weight was observed in the HAXPES. In addition, $c$-axis elongation also increased spectral weight away from the Fermi energy (i.e., at $∼$1.5 eV). Our simulations presented in Fig. 6 suggest that the LHB contributes to this region and is dependent on the strain orientation. In order to interrogate these features further, we have employed resonant x-ray photoemission spectroscopy (RPES) at the V L$3$-edge to enhance our sensitivity of the occupied V 3$d$ orbital region. The 10 nm epitaxial films were grown by reactive molecular beam epitaxy at Cornell University, further details are provided by Paik et al.¹² Fresh samples were prepared ahead of the RPES experiments at the 29-ID beamline at Argonne Photon Source on the soft x-ray ARPES end-station. The samples were measured in their respective high temperature phases [at room temperature for (001) and at 398K for (100)], and spectra were energy calibrated using a combination of Au foil and VO$2$ film references. When at resonance, the direct recombination of the decay process emits 3$d$ electrons when excited at the L-edge, resulting in enhanced 3$d$ contributions in the measured x-ray photoemission spectra. Eguchi et al. previously employed V L-edge RPES to measure the insulating and metallic phases of 10 nm thin VO$2$ (001) films, where they were able to enhance sensitivity to the quasiparticle peak and incoherent feature.⁴² For the VO$2$ (001) case, we reproduced the metallic phase RPES spectra presented in Ref. 42 for the same energies they employed, as shown in Fig. 7(b). Figure 7(c) directly compares the resonant spectra for the VO$2$ (001) and (100) cases for select energies. The largest enhancement compared to the off-resonance data was observed for the energies resonant with the knee and main peak of the V L$3$-edge absorption line. We note that these energies displayed the largest strain-induced dichroism due to the preferential orbital filling associated with the strain-induced OSMT.¹⁴ In addition to the weaker quasiparticle peak intensity of the VO$2$ (100) compared to the (001) case that we reported previously,¹⁴ the increased sensitivity afforded by the resonance reveals how this effect is correlated with increased incoherent peak at $∼$1.5 eV for the (100) case. These data agree with simulated spectral functions in Fig. 6, where a stronger incoherent peak at $∼$1.5 eV is always observed for the case where the c-axis is elongated, i.e., (100), compared to bulk or when the c-axis is compressed, i.e., (001). As a result, we can assign this feature as the lower Hubbard band and conclude that VO$2$ can be strain-engineered as predicted by our modeling of the OSMT. Figure 8 summarizes the general trends of the redistribution of electron occupation numbers due to the interplay between strain and Hubbard interactions in the epitaxial VO$2$/TiO$2$ systems.

As a comparison, in this section, we study the strain effect on NbO$2$ in the rutile phase. NbO$2$ undergoes a metal-to-insulator transition at a temperature around 810 $°$C. It simultaneously undergoes a crystal structure change between a high temperature, undistorted rutile phase and a low temperature, body-centered tetragonal, distorted rutile phase.⁴³ Although strain engineering on NbO$2$ in the rutile phase has not yet been experimentally achieved, here we investigate the change of quasiparticle weights in rutile NbO$2$ under the same strain condition as the VO$2$ case. For the bulk NbO$2$⁠, we adopt the LDA-optimized lattice parameters, $(aR,bR,cR)=(4.93$Å,4.93Å,2.9Å).⁴³ For NbO$2$ (100), we choose the lattice parameters of $(aR′,bR′,cR′)=(4.74$Å,4.93Å,3.016Å), which is in the same strain condition of VO$2$ (100) studied in Sec. V. Note that the spin-orbit coupling is not included in our calculations. Figures 9(a) and 9(b) plot the non-interacting DOS obtained from the DFT-fitted $HTB$ for bulk and (100), respectively, and the orbital-dependent quasiparticle weights calculated from our DFT-SS method with $U=6$ eV and $J=1$ are summarized in Figs. 9(c) and 9(d). We find that although the bandwith of the $dx2−y2$ orbital still has a significant reduction in (100) compared to the bulk, the quasiparticle weights are almost the same in NbO$2$ with and without the strain, which is very different from the VO$2$ case. We attribute this result to the fact that the bandwidth of $d$ orbitals in NbO$2$ is larger than that of VO$2$⁠, and thus, the ratio of the interaction energy to the bandwidth (⁠$U/W$⁠) is generally smaller even with the same values of interaction parameters. This suggests that NbO$2$ has a weaker correlation and, consequently, is not an ideal material for modulating the Mott physics with strain.

To experimentally investigate this conclusion, HAXPES was performed on rutile VO$2$ (100), VO$2$ (001), and relaxed NbO$2$ (110) thin films, as shown in Fig. 10. The HAXPES measurements were performed with 5.95 keV x-rays at the I09 beamline of the Diamond Light Source. HAXPES spectra were energy-resolved and measured using a Scienta Omicron EW4000 high-energy analyzer with a $±30°$ acceptance angle. The x-ray beam was monochromatized using a Si (111) double-crystal monochromator, followed by a Si(004) channel-cut crystal resulting in an overall energy resolution of $≤250$ meV. The films were all heated above their respective MIT temperatures during measurement to transition them into their metallic, rutile phases. It can be clearly seen that unlike the VO$2$ thin films, the Nb $d$-orbital feature near the Fermi level of NbO$2$ can be well-fit by a single peak without the need for an additional hump feature associated with the lower Hubbard band. This observation is completely consistent with our conclusion on the NbO$2$ discussed above.

In this paper, we have studied general trends in the redistribution of electron occupation numbers due to the interplay between strain and correlations in epitaxial VO$2$/TiO$2$ systems, which is summarized in Fig. 8. Since the $eg$ bands are far from the Fermi energy, the correlation effects will be related mainly to the $t2g$ bands. We find that a stronger Hubbard interaction tends to make the electron fillings in $t2g$ orbitals more even, which is a general trend for a non-zero Hund’s coupling $J$⁠. Moreover, because most extra electrons will go into the $t2g$ orbitals as $ntot$ increases, increasing $ntot$ will push the system closer to the half-filling condition (⁠$n=3$ for a three-orbital $t2g$ system). This consequently enhances correlation effects, resulting in the reduction of quasiparticle weights, and the emergence of the incoherent spectrum as shown in our calculations. Finally, the strain in VO$2$ (100) causes an elongation of the rutile $c$-axis that produces a significant reduction of the bandwidth, particularly on the $d∥$ band, but not on the $dπ$ bands. This orbital-dependent reduction of the bandwidth leads to an intriguing “orbital-selective” Mott transition (OSMT) in which only one of the orbitals is pushed much closer to the Mott insulating state. Based on our results, we conclude that in epitaxial VO$2$/TiO$2$ (100) and (110),¹⁴ the system moves toward the OSMT due to $c$-axis elongation by the strain from lattice matching with the TiO$2$ substrate, as indicated by the green arrow in Fig. 8.

In summary, we have employed density functional theory combined with the slave spin method (DFT+SS) to study the electronic structures of epitaxial VO$2$ films under strain in the rutile phase. We have found that while asymmetrical orbital occupation numbers are favored due to band structure effects in the strained systems with an elongated $c$-axis, strong electron-electron correlation drives orbital-dependent modifications of the quasiparticle weights, as well as the incoherent spectra in the single particle Green’s function. The interplay between these two distinct effects pushes epitaxial VO$2$ films toward an orbital selective Mott transition (OSMT) in the rutile phase without a Peierls instability. Our results indicate that the Mott physics is important and can be significantly modulated even in the rutile phase of VO$2$ by appropriate strain-engineering.

This work is supported by the Air Force Office of Scientific Research under Award No. FA9550-18-10024. This research used the resources of the Advanced Photon Source, a U.S. Department of Energy (DOE) Office of Science User Facility operated by the DOE Office of Science by Argonne National Laboratory under Contract No. DE-AC02-06CH11357; additional support was provided by the National Science Foundation under Grant No. DMR-0703406. Part of the calculations used the Extreme Science and Engineering Discovery Environment (XSEDE) supported by the National Science Foundation (Grant No. ACI-1548562). A.R. acknowledges the support from Imperial College London for her Imperial College Research Fellowship. We thank Diamond Light Source for access to beamline I09 (SI20647 and SI21430-1) that contributed to the results presented here. J.E.B. acknowledges funding by the Leibniz association within the Leibniz Competition through a project entitled “Physics and control of defects in oxide films for adaptive electronics.”

The motivation behind the slave-spin formalism is to introduce a two-level degree of freedom to handle the charge dynamics of interacting fermions. The occupation of each local degree of freedom, whether it be site, spin, or orbital, can be treated as either occupied (1) or unoccupied (0), so a spin algebra can be used to represent this. That means, the spin operator seen in the following derivations is not for the physical spin, but rather, it is a “slave” spin capturing the dynamics of “charge” degrees of freedom. The physical spin is described by fermionic operators called spinons. With these definitions, one can rewrite the electron creation operator as

where $i$ is the site index, $α$ is the orbital index, and $σ$ is the physical spin. The spinon $fiασ†$ is a fermionic operator capturing the dynamics of the physical spin. Because the Hilbert space for the charge degree of freedom is parameterized as two levels (states with 0 and $−e$ charge), which is fully described by a spin 1/2 algebra, we can have $Siασ+$ as the spin raising operator for the charge part of the electrons. Since we have a slave spin for each set of quantum numbers $(α,σ)$⁠, $Siασ+$ carries both the orbital and the physical spin indices. Similar to other slave particle formalisms, the Hilbert space has been enlarged in this process and will contain unphysical states. For example, a state such as $|niασf=0,Siασz=+1/2⟩$ would be unphysical and must be projected out of the enlarged space. This projection can be done by enforcing the constraint

Following the procedure outlined in Ref. 31, we can define the mean-field parameters

and $λα,σ$ being the Lagrange multipliers to enforce the constraint on the average over the sites

where $N$ is the total number of sites. At the mean-field level, we have the decoupled Hamiltonians for the spinons (⁠$Hf$⁠) and for the slave spins (⁠$HS$⁠), which can be used to solve ${zα,σ,λα,σ}$ self-consistently. After the self-consistent solution is obtained, the quasiparticle weight for each orbital can be evaluated by

Because we are only interested in the paramagnetic phase, we have $Zα=|zα,↑|2=|zα,↓|2$⁠.

Using the same procedure given in Ref. 30, the imaginary-time ordered Green function can be evaluated at the mean-field level as

In the Lehmann representation, the retarded Green function can be evaluated directly by

where ${|m⟩}$ and ${Em}$ are the sets of eigenvectors and eigenvalues of the slave spin mean-field Hamiltonian $HS$⁠, $λ$ is the band index, $nλf(k→)$ is the Fermi distribution function for the state with momentum $k→$ in $λ$th band, ${ϵλ(k→)}$ is the set of eigenvalues of the spinon mean-field Hamintonian $Hf$⁠, and $Uk→σij$ is the unitary matrix that diagonalizes $Hf$⁠. Equation (B2) can be evaluated directly without any difficulty after the mean-field equations are solved.

With the electron Green function given in Eq. (B2), the electron spectral function for electrons in each orbital can be obtained straightforwardly by

One common problem with the slave particle approaches is the missing spectral weight at doping away from the half-filling. We check the sum rule

and we find that the integrated spectral weight equals to 1 only as $2nασ−1=0$⁠, which corresponds to the “half-filling” in the orbital $α$⁠. This failure is an artifact of the mean-field approximation in which the couplings between spinons and slave spins are neglected in Eq. (B1). It is generally expected that once the couplings are turned on, the spinon Fermi surface will be smeared out. As a result, we can set $nσk→f=1/2$ in Eqs. (B2) and (B4) at the mean-field level, which leads to the correct sum rule of

Finally, as pointed out in Ref. 30, the contribution to the spectral function from $m=n$ in Eq. (B2) is the coherent part, while other contributions from $m≠n$ correspond to the incoherent part.