Charge disproportionation and Hund’s insulating behavior in a five-orbital Hubbard model applicable to $d^{4}$ perovskites

We study the following model, which consists of a single-particle tight-binding (TB) part and a local interaction on each site $R$:

The TB part consists of inter-site hopping terms as well as an on-site crystal-field splitting, $\Delta_{{e_{g}}-{t_{2g}}}$, between ${e_{g}}$ and ${t_{2g}}$ orbitals. We include hopping between nearest neighbors (NN) and next-nearest neighbors (NNN) on a cubic lattice, representing the sites of the transition-metal cations in the perovskite structure. The hopping amplitudes are assumed to have the cubic symmetry corresponding to ${t_{2g}}$ and ${e_{g}}$ orbitals. Thus, between nearest neighbors there is only hopping between ${t_{2g}}$ orbitals of the same type and no hopping connecting the ${e_{g}}$ and ${t_{2g}}$ orbitals. For simplicity, we neglect next-nearest neighbor hopping between different $t_{2g}$ and between $t_{2g}$ and $e_{g}$ orbitals, and we assume the same fixed ratio $t_{\mathrm{NNN}}/t_{\mathrm{NN}}$ between the NNN and NN hopping amplitudes for both ${e_{g}}$ and ${t_{2g}}$ orbitals. For the $e_{g}$ orbitals, $t_{\text{NN}}$ and $t_{\text{NNN}}$ correspond to the hopping between $d_{3z^{2}-r^{2}}$ orbitals along $z$ and along $x+z$ or $y+z$, respectively, see Ref. 33.

We also perform calculations where we introduce an intrinsic “site splitting” of the $e_{g}$ levels, by raising, respectively lowering, the local on-site energy of the ${e_{g}}$ orbitals on adjacent sites by $\pm\Delta_{e_{g}}/2$. This mimics the effect of the structural $R_{1}^{+}$ breathing mode distortion that typically accompanies the charge disproportionation in perovskites (see, e.g., Ref. 34). As described in Sec. I, this distortion consists in a spatially alternating expansion and compression of the oxygen octahedra surrounding the transition-metal sites, resulting in a three-dimensional checkerboard-like pattern of large and small octahedra.

To describe the local electron-electron interaction within the five-orbital $d$ shell, we use the so-called Slater parametrization in the density-density approximation (see, e.g., Ref. 35):

This term is purely local and identical for each transition-metal site (we therefore suppress the site index $R$ in Eq. (II.1)). Thereby, $n_{m\sigma}$ is the occupation number operator for orbital $m$ with spin $\sigma$ (on site $R$) and $\bar{\sigma}$ is short for $-\sigma$. The matrix elements $U_{mm^{\prime}}$ and $J_{mm^{\prime}}$ are parametrized in terms of the Slater integrals $F_{0}$, $F_{2}$, and $F_{4}$ in the usual form [35]. The Slater integrals themselves are expressed in terms of the average Coulomb repulsion $U=F_{0}$ and Hund’s rule interaction $J=(F_{2}+F_{4})/14$, using a fixed ratio $F_{4}/F_{2}=0.63$ [35].

We also compare the full five-orbital $d$-shell model to a simplified effective two-orbital ${e_{g}}$ model, applicable to the high-spin limit. In this model, the occupation of the ${t_{2g}}$ states on each site is fixed to three electrons, giving rise to an $S=3/2$ core spin, and hopping in and out of these ${t_{2g}}$ states is neglected [31; 32]. The remaining electron is thus constrained to the $e_{g}$ orbitals, but can still hop from site to site. Then, for every site, the ${t_{2g}}$ core spin defines a local spin quantization axis, and the spins of the ${e_{g}}$ electrons are defined relative to this local axis, denoted by $\sigma=\Uparrow$/$\Downarrow$.

The $e_{g}$ electrons interact with the ${t_{2g}}$ core spins through an effective Hund’s rule interaction, which leads to a local spin splitting of the ${e_{g}}$ orbitals:

where $h$ is the effective magnetic field created by the ${t_{2g}}$ core spin. The strength of this effective magnetic field can be related to the parameters of the full interaction Hamiltonian, Eq. (II.1), in particular to the Hund’s interaction parameter $J$, by considering the energy gain (penalty) of an ${e_{g}}$ spin (anti-)parallel to the ${t_{2g}}$ core spin, i.e., by considering the energy difference between $(t_{2g})^{3}_{\uparrow}(e_{g})^{1}_{\uparrow}$ and $(t_{2g})^{3}_{\uparrow}(e_{g})^{1}_{\downarrow}$:

Furthermore, assuming a completely random, i.e., paramagnetic orientation of the ${t_{2g}}$ core spins on the different sites results in a renormalization of the hopping amplitudes, enabling also hopping between $\Uparrow$ and $\Downarrow$ states corresponding to different sites. After averaging over all possible relative spin orientations, this can be described by a matrix $u^{R,R^{\prime}}_{\sigma,\sigma^{\prime}}$, which does not affect the on-site Hamiltonian ($u^{R,R}_{\sigma,\sigma^{\prime}}=\delta_{\sigma,\sigma^{\prime}}$) but renormalizes and mixes the inter-site hopping terms ($u^{R,R^{\prime}\neq R}_{\sigma,\sigma^{\prime}}=2/3$) according to [31; 32]:

To enable a consistent comparison with the full five-orbital model, the interaction in the effective two-orbital model is described using the ${e_{g}}$ subspace of the full interaction Hamiltonian, Eq. (II.1). This is equivalent to the density-density approximation of the so-called Kanamori Hamiltonian [3], but with the corresponding interaction parameters defined in terms of the Slater integrals (or in terms of the average Coulomb repulsion $U$ and the corresponding Hund’s rule interaction $J$).

The purpose of this work is to study the physics of charge disproportionation in a relatively generic cubic five-orbital model. Nevertheless, we focus on realistic parameter regimes, applicable to perovskite transition-metal oxides and closely related materials. We therefore obtain parameters for the single-particle part of the Hamiltonian, $H_{\text{TB}}$, by comparison with the full band structure of CaFeO₃, calculated within density functional theory (DFT).

We perform DFT calculations using the “Vienna Ab-initio Simulation Package” (VASP) [36; 37] employing the PBEsol exchange-correlation functional [38], both for the high-temperature $Pbnm$ structure, where all Fe sites are equivalent, and the charge-disproportionated $P2_{1}/n$ crystal structure observed below $290\text{\,}\mathrm{K}$ [39; 27]. We use a $7\times 7\times 5$ $k$-point mesh for the 20-atom unit cell and a plane-wave energy cutoff of $600\text{\,}\mathrm{e}\mathrm{V}$ to ensure well converged results.

We first relax unit-cell parameters as well as atomic positions within $Pbnm$ symmetry. In order to stabilize the high-spin state of the Fe cation, a Hubbard-corrected DFT+$U$ treatment and the presence of antiferromagnetic (AFM) order is required. We find that moderate values of $U=$4\text{\,}\mathrm{e}\mathrm{V}$$ and $J=$1\text{\,}\mathrm{e}\mathrm{V}$$ and A-type antiferromagnetic order (as proxy for the more complicated helical spin structure in CaFeO₃ [40]) result in structural parameters that are in good agreement with experimental data (e.g., the calculated unit cell volume is $210.34\text{\,}{\mathrm{\SIUnitSymbolAngstrom}}^{3}$) [39; 26]. We then re-relax the atomic positions within non-spin-polarized DFT (to exclude any structural effects or symmetry lowering stemming from the magnetic order) but with lattice parameters fixed to the ones obtained from the magnetic $+U$ calculations. We note that a non-magnetic DFT calculation leads to a low-spin state of the Fe cation, which would result in an unrealistically low unit-cell volume.

To obtain a realistic estimate of the site-splitting $\Delta_{e_{g}}$, we manually add a breathing mode of ${R_{1}^{+}}=$0.2\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$, very similar to the experimentally observed amplitude ${R_{1}^{+}}=$0.18\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$ [26], while keeping all other structural parameters unchanged. A full DFT+$U$ study of CaFeO₃ has been presented in Ref. 41.

The calculated band structure and corresponding densities of states (DOS) in the relevant energy window around the Fermi level are shown in Fig. 1. It can be seen that a group of bands with dominant Fe character, located between 0.5 eV and 4 eV above the Fermi level, is separated from all other bands at higher and lower energies. Further inspection confirms that these are the Fe-${e_{g}}$ bands. The bands between approximately $-$1 eV and 0.5 eV also have mostly Fe character and correspond to the Fe-${t_{2g}}$ states. Towards lower energies, there is some entanglement of these Fe-${t_{2g}}$ bands with other bands with dominant O-$p$ character. Comparing the band structures obtained with and without the structural $R_{1}^{+}$ breathing mode, one can see that its main effect is to introduce a splitting throughout the middle of the ${e_{g}}$ bands, analogously to what has been observed for nickelates [8; 42]. In contrast, the Fe-${t_{2g}}$ bands are not much affected by the structural distortion.

We fit the parameters of our TB model such that the corresponding band dispersion matches the bands obtained from our DFT calculations. Thereby, the ratio between NNN and NN hopping, $t_{\text{NNN}}/t_{\text{NN}}$, is obtained by fitting the well-separated ${e_{g}}$ bands, whereas the ${t_{2g}}$ NN hopping is then mainly adjusted to obtain the correct ${t_{2g}}$ band width. Note that the two cases with and without breathing mode are fitted with the same hopping amplitudes, whereas the average ${e_{g}}$-${t_{2g}}$ splitting $\Delta_{{e_{g}}-{t_{2g}}}$ is slightly altered by the ${R_{1}^{+}}$ breathing mode. The resulting TB band structure is also shown in Fig. 1, and the corresponding parameters are listed in Table 1. It can be seen that the TB model results in an excellent fit of the Fe-${e_{g}}$ bands, both in the $Pbnm$ and $P2_{1}/n$ cases. For the ${t_{2g}}$ bands, the overall band width is well described, while the dispersion of the individual bands shows notable deviations. This could be improved by introducing a separate $t_{\text{NNN}}/t_{\text{NN}}$ for the ${t_{2g}}$ hopping or other additional intra- and inter-site TB parameters. However, for simplicity, and since the specific form of the ${t_{2g}}$ bands is not relevant in the context of this work, we refrain from introducing further parameters into our TB model.

We obtain the phase diagram of the model, Eq. (1), as function of the interaction parameters $U$ and $J$ using dynamical mean-field theory (DMFT) [43]. Within DMFT, the self-energy of the full lattice problem is approximated by a frequency-dependent but local self-energy, which is obtained by mapping each symmetry-inequivalent site on an effective impurity problem, defined via a self-consistency condition requiring that the impurity Green’s function is identical to the corresponding local Green’s function of the full lattice problem.

The calculations are implemented using the TRIQS and DFTTools libraries [44; 45; 46]. To solve the effective impurity problem, we use a continuous-time Quantum-Monte-Carlo solver (CT-HYB), also implemented in the TRIQS library [47], and an inverse temperature $\beta=1/(k_{\text{B}}T)=$40\text{\,}\mathrm{e}\mathrm{V}^{-1}$$, corresponding to approximately room temperature. Green’s functions are represented by 40 Legendre coefficients [48] for the full five-orbital model (50 for the effective two-orbital model), which results in good convergence and a smooth behavior of the self-energy in imaginary frequency $\Sigma(\mathrm{i}\omega_{n})$.

We average over both spin channels and also over the ${e_{g}}$ and ${t_{2g}}$ manifolds separately to ensure a paramagnetic solution without orbital polarization in the two different subshells. Furthermore, to allow for a charge disproportionated solution, we divide the sites of the simple cubic lattice into two interpenetrating fcc lattices, (I) and (II), corresponding to the expected 3D checkerboard-like arrangement, and solve the two resulting effective impurity problems separately.

From the local Green’s function $G_{mm^{\prime}}(\tau)=-\langle Tc_{m}(\tau)c^{\dagger}_{m^{\prime}}(0)\rangle$, where $T$ is the imaginary-time ordering operator and $c^{\dagger}_{m}(\tau)$ is the creation operator for an electron in orbital $m$ at the local site at imaginary time $\tau$, we obtain the total ${e_{g}}$ occupation per site, $n_{\mathrm{e_{g}}}=n_{\mathrm{z^{2}}}+n_{\mathrm{x^{2}-y^{2}}}$, the (positively-defined) occupation difference between the two sub-lattices, $\delta n=|\sum_{m}n_{m}^{\mathrm{(I)}}-\sum_{m}n_{m}^{\mathrm{(II)}}|$, and the averaged spectral weight around the Fermi level, $A(\omega=0)=-(\beta/\pi)\mathrm{Tr}\,G(\tau=\beta/2)$. Furthermore, we obtain an estimate for the (orbital-dependent) quasiparticle weight $Z_{m}=\left[1-\left.\frac{\partial\imaginary\Sigma_{m}(\mathrm{i}\omega)}{\partial(\mathrm{i}\omega)}\right|_{\mathrm{i}\omega=0}\right]^{-1}$ by fitting a fourth-order polynomial to the imaginary part of the local self-energy, $\imaginary\Sigma_{m}(\mathrm{i}\omega_{n})$, for the eight lowest positive Matsubara frequencies $\mathrm{i}\omega_{n}$ and extrapolate the slope at zero frequency from this fit, c.f. Appendix C in Ref. 49. We use the maximum-entropy method [50], implemented in the TRIQS library [51], to obtain both $k$-averaged and $k$-resolved spectral functions on the real frequency axis.

Before presenting our results from the full DMFT calculations, we outline the general features of the expected phase diagram by considering the “atomic limit” of Eq. (1), i.e., the limit where all hopping amplitudes and thus the associated band widths go to zero, and the ground state configuration of the system is solely determined by the local part of the Hamiltonian, i.e., the interaction term $H_{\text{int}}$ plus the local crystal-field and ${e_{g}}$ site splitting.

We first compare the energies of the three most relevant configurations: (i) homogeneous high-spin (HS) $d^{4}$; (ii) homogeneous low-spin (LS) $d^{4}$; and (iii) inhomogeneous HS $d^{5}$ + $d^{3}$, also depicted by the schematic level diagrams in Fig. 2, and identify the ground state as function of interaction parameters $U$ and $J$. Details are described in App. A. Fig. 2 depicts the case without an explicit ${e_{g}}$ site splitting, i.e., $\Delta_{e_{g}}=0$. Thick black lines separate the regions where the three different cases (i)-(iii) are the ground state configuration. The homogeneous LS configuration (ii) is generally favorable for small $J$ ($J\lesssim\Delta_{e_{g}-t_{2g}}/2.83\approx 0.88$ eV, except for the small part separating the LS and inhomogeneous HS regions, where the critical $J$ is further decreased by decreasing $U$). The inhomogeneous charge-disproportionated (CD) state (iii) is the ground state for larger $J$ and $U\lesssim 1.51J$ (except, again, for the small part separating CD and homogeneous LS regions, where the critical $U$ is smaller). For large $J$ and $U$, the homogeneous HS state (i) is obtained. We note that in the Slater parametrization of the local interaction, a positive $U$ represents an average repulsive electron-electron interaction (independent of $J$), which is why we consider the whole region with $U>0$ as potentially physically relevant.

A simple estimate for the boundary between metallic and (Mott-) insulating regions can be obtained by comparing the lowest inter-site charge transfer excitation energy with the average electronic kinetic energy. For the homogeneous phases, this kinetic energy can be approximated by the relevant band width (see, e.g., Ref. 4), resulting in the blue-shaded “Mott-insulating” region in Fig. 2 (see App. A for details). Since the ${t_{2g}}$ band width is smaller than the ${e_{g}}$ band width, the Mott phase extends to smaller $U$ values in the LS region than in the HS region.

For the CD region, a good estimate for the corresponding effective band width/kinetic energy is more difficult to obtain. For the transition to the homogeneous metallic phase it turns out to be significantly smaller than $W_{e_{g}}$ and essentially goes to zero when $\Delta_{e_{g}}$ becomes large [8]. For simplicity, we therefore denote the whole region below the CD transition line as “CDI” in the atomic limit (red shaded area in Fig. 2). Considering only hopping among either nominal $d^{5}$ or $d^{3}$ sites defines limits for the CDI region towards small $U$ (see App. A). Simple band-width estimates for the kinetic energy result in the limits indicated by the red-striped regions in Fig. 2.

Introducing a nonzero ${e_{g}}$ site splitting $\Delta_{e_{g}}$ (corresponding to a nonzero ${R_{1}^{+}}$ distortion of the underlying crystal structure), extends the CD region slightly by rigidly shifting the phase boundaries towards the homogeneous HS and LS phases (see App. A) as indicated by the small arrows in Fig. 2. This also shifts the metal/insulator phase boundary within the homogeneous HS region towards higher $U$. In contrast, reducing the ${e_{g}}$ band width shifts the metal/insulator phase boundary in the homogeneous HS phase in the opposite direction, while not affecting any of the other main phase boundaries.

Next, we present the results of our full DMFT calculations for the TB model with $\Delta_{e_{g}}=0$, i.e., for the case without a site-dependent shift of the $e_{g}$ levels. Fig. 3 depicts various observables as function of interaction parameters $U$ and $J$.

The total $e_{g}$ occupation per site, $n_{e_{g}}$, shown in Fig. 3(a), indicates a sharp transition from LS ($n_{e_{g}}=0$) to HS ($n_{e_{g}}=1$) for $J\geq 1$ eV, essentially independent of $U$, in very good agreement with the expectation obtained from the atomic limit discussed in Sec. III.1. Only for $J=0.75$ eV and $U\leq 2$ eV, we obtain a small intermediate-spin region where the $e_{g}$ occupation lies between 0 and 1.

The most important result is visible in the occupation difference between the two sublattices, $\delta n$, shown in Fig. 3(b). Here, in the region of small $U$ but large $J$, i.e., in the region where the atomic limit predicts the inhomogeneous CD ground state, we indeed obtain a non-zero value for $\delta n$, indicating that the system exhibits a spontaneously symmetry-broken state with higher and lower $e_{g}$ occupation on the two sublattices. Note that this CD region lies fully within the HS region where the average $e_{g}$ occupation is equal to 1. Inspection of the spectral weight at zero energy, $A(0)$, shows that the charge disproportionation coincides with a transition to an insulating state, indicated by $A(0)\approx 0$ (dark regions in Fig. 3(c)). This CD insulating region is separated from the Mott-insulating regime, where $\delta n=0$ and which occurs for large $U$, by a metallic region. It can be seen that the boundary between the Mott-insulating and the metallic region (both in the HS and LS regime) also matches very well with the atomic limit prediction from Sec. III.1, indicated by the gray lines in Fig. 3. On the CD side, the MIT occurs together with the emerging charge disproportionation for $\delta n\gtrsim 1$, somewhat below the homogeneous-inhomogeneous transition line predicted from the atomic limit. However, we note that the distance of the atomic-limit prediction and the actual charge disproportionation is significantly less than the $e_{g}$ band width.

In Fig. 3(c), there is a significant change in spectral weight, $A(0)$, at the LS/HS transition in the metallic region. Fig. 4(a) shows corresponding spectral functions (for $U=1.5$ eV), separated in $t_{2g}$ and $e_{g}$ contributions. It can be seen that for $J=0.5$ eV the full spectral weight is generated by the $t_{2g}$ contribution, while the $e_{g}$ states are still unoccupied. Around $J=0.75$ eV the transition occurs, where both $e_{g}$ and $t_{2g}$ states contribute to the spectral weight at zero energy. Above this LS/HS transition, at $J=1$ eV, the $t_{2g}$ spectral function is gapped, and $A(0)$ is entirely due to the $e_{g}$ states, i.e., a situation that corresponds to an orbitally selective Mott state. Here, the $e_{g}$ quasiparticle peak is lower and broader compared to the $t_{2g}$ peak for $J=0.5$ eV, which results in the lower spectral weight in the phase diagram for the HS case.

Fig. 3(d) shows the quasiparticle weight $Z$ in the metallic region, averaged over all orbitals with non-vanishing spectral weight. $Z$ is largest along the line separating homogeneous/inhomogeneous ground states in the atomic limit, and decreases towards both insulating phases. This is consistent with analogous results for the three-orbital case in Ref. 9, where this behavior was related to the competition between the two different insulating phases, resulting in the so-called “Janus effect” [52; 4].

Fig. 4(b) depicts the imaginary part of the self-energy on the Matsubara axis across the LS/HS transition within the metallic regime. As described in Sec. II.4, the quasiparticle weight $Z$ is obtained by fitting the low-frequency behavior of $\imaginary\Sigma(\mathrm{i}\omega)$. The corresponding fits are shown in Fig. 4(b). It can be seen that the self-energies exhibit rather strong deviations from simple Fermi liquid behavior. A linear dependence on $\mathrm{i}\omega$ can only be observed in a very small range close to $\omega\rightarrow 0$, and there is a rather large non-vanishing imaginary part for zero frequency. This could in part be due to the relatively high temperature ($\beta=40\,$eV, corresponding to approximately room temperature) in our calculations. Thus, the interpretation of $[1-\partial\imaginary\Sigma/\partial(\mathrm{i}\omega)]^{-1}$ as effective quasiparticle renormalization is somewhat approximate, here. This can also be inferred from the comparison with the “quasiparticle weights” extracted by integrating over the central peak of the corresponding spectral functions in Fig. 3. (Note that the latter of course also just represents a very rough and approximate measure of the quasiparticle renormalization.)

To obtain further insight into the different phases, we also analyze the probabilities of the different atomic multiplet states, obtained from the DMFT Quantum-Monte-Carlo solver. These probabilities describe how likely it is for the corresponding state to be visited during the simulation of the effective impurity problem and allow to distinguish between different occupations and LS/HS configurations in a more quantitative way than in the phase diagrams shown in Fig. 3. To simplify the analysis, we sum up the probabilities for Hilbert states with the same number of electrons $N$ and the same magnetic spin quantum number $|m_{s}|$ ¹¹1We note that within the density-density approximation for the interaction Hamiltonian, Eq. (II.1), the magnetic quantum number $m_{s}$ and not the total spin $S$ is the relevant quantum number indicating the spin state of the system.. It turns out that only states with $N=3$, $4$, and $5$ have significant probabilities; these together with the quasiparticle weight $Z$ are shown along a rectangular path through the phase diagram in Fig. 5.

The first segment of this path corresponds to constant $U=1.5$ eV and increasing $J$ from 0.5 eV to 2.0 eV, i.e., starting from the LS metallic, then to the HS metallic, and then into the CDI phase. In the LS metallic phase for $J=0.5$ eV, all multiplets with $N=3$, $4$, and $5$ are visited, but the probabilities of the HS states (maximal $|m_{s}|$) for $N=4$ and 5 are close to zero. In the transition region between LS and HS metallic phases at $J=0.75$ eV, the HS multiplets with $|m_{\mathrm{s}}|=N/2$ for $N=4$ and $N=5$ become more populated, while for $J=1$ eV, the HS states amount to nearly the total probabilities.

For increasing $J$, the probability for $N=4$ decreases linearly while the probabilities for $N=3$ and $N=5$ increase accordingly. When the system enters the CDI phase ($J>1.5$ eV), the probabilities on the two inequivalent sites become different, but the averages over both sites, shown in Fig. 5, continue their linear trend. In particular, even deep inside the CDI phase, for $J=2$ eV, the probability of $N=4$ remains relatively large on both sites, indicating significant fluctuations between multiplets with different $N$, even in the insulating state. We note that the presence of such fluctuations and the finite probability of the $N=4$ multiplet is consistent with the gradual increase of the occupation difference $\delta n$ within the CDI phase shown in Fig. 3.

Fig. 6 shows the site-resolved ${e_{g}}$ spectral functions in this regime ($U=1.5$ eV and $J\geq 1.5$ eV). The spectral functions for $J>1.5$ eV exhibit a clear gap, consistent with the zero spectral weight in Fig. 3, despite the system not being in a pure local charge state, as follows from the multiplet probabilities. This differs from a recent study using a slave-boson method [9], where the Hund’s insulator in the two-orbital model was characterized by a complete charge disproportionation where only the nominal charge multiplets are populated. It also indicates that the transition to the CDI state cannot be understood purely in terms of atomic-limit considerations such as the ones in Sec. III.1.

Increasing $U$ from 1.5 eV to 6.0 eV for fixed $J=1.5$ eV (second segment in Fig. 5) results in a transition from the CDI phase back to the HS metallic region and then moves the system towards the Mott insulating region. Here, the $N=4$ probability increases again, while the $N=3$ and $N=5$ probabilities decrease accordingly, and become zero at the boundary to the Mott insulator. Also, the quasiparticle weight $Z$ decreases continuously towards zero when approaching the Mott phase, while it remains quite large even very close to the CDI phase boundary. However, as discussed above, it is unclear whether $[1-\partial\imaginary\Sigma/\partial(\mathrm{i}\omega)]^{-1}$ is really a good quantitative measure for the quasiparticle renormalization, here. Nevertheless, we note that a rather sharp drop of $Z$ on approaching the Hund’s insulator (equivalent to our CDI phase) has also been observed in Ref. 9 for a three-orbital model.

Finally, when crossing the HS/LS transition within the Mott insulating region, the HS $N=4$ multiplet vanishes abruptly, and then, when leaving the Mott phase towards the LS metallic region, the $N=3$ and $4$ multiplets appear again, together with a gradual increase of $Z$.

Thus, to preliminary summarize the results in this section, we emphasize that we establish the presence of a spontaneously charge-disproportionated insulating phase, for strong $J$ and intermediate values of $U$, in the cubic five-orbital multi-band Hubbard model with an average of four electrons per site, analogous to results obtained for similar two- and three-orbital models [8; 9]. The CDI phase is characterized by a gap opening at zero frequency, but still exhibits a strong local mixing of multiplets with different number of electrons $N$. This is in contrast to the Mott insulating phase, which can be characterized by a unique local multiplet with $N=4$. In addition to the metallic and two distinct insulating phases, a LS/HS transition occurs. This transition is very sharp inside the Mott-insulating regime but a small intermediate spin region can be observed for small $U$ in the metallic phase.

Next, we investigate the effect of an intrinsic energy difference between the ${e_{g}}$ levels on the two sublattices, which is introduced, for example, by the structural distortion that generally accompanies the charge disproportionation, i.e., the ${R_{1}^{+}}$ breathing mode of the oxygen octahedra surrounding the transition-metal sites in the perovskite structure. We use a splitting of $\Delta_{e_{g}}=0.436$ eV as obtained from our fit to the band structure of CaFeO₃ described in Sec. II.3.

Fig. 7 shows the spectral weight $A(0)$ and occupation difference $\delta n$ for this case. It can be seen that the region corresponding to the CDI phase is significantly increased compared to the case with $\Delta_{e_{g}}=0$ in Fig. 3. The transition to the insulating state now occurs already very close to the homogeneous/inhomogeneous transition line estimated from the atomic limit. Similar behavior has also been observed for the two-band model [8]. Furthermore, due to the intrinsic energy difference, a small $\delta n>0$ already develops in the metallic phase. On the other hand, this small occupation difference is completely suppressed when approaching the Mott-insulating phase, and the corresponding MIT occurs again very close to the atomic limit predictions based on the $e_{g}$ band width. Note that the corresponding lines in Fig. 7 are shifted relative to those in Fig. 3 for $\Delta_{e_{g}}=0$ as indicated in the schematic phase diagram, Fig. 2.

Due to the large extent of the CDI phase, there is now a direct competition between the CDI and the homogeneous LS metallic phase for small $J$. Nevertheless, also in this case, the corresponding phase boundary appears to be rather well described by the atomic limit prediction (gray lines in Fig. 7). Similar to the case with $\Delta_{e_{g}}=0$, we obtain a small intermediate-spin region (without significant charge disproportionation) for $U\leq 2$ eV and $J=0.75$ eV (data not shown in Fig. 7).

Finally, we also assess the effect of reducing the $e_{g}$ band width, while keeping the local crystal field and the inter-site ${e_{g}}$ splitting constant. To simplify the corresponding analysis, we use the observation that within the HS region ($J\geq 1$ eV), the gradient of all quantities depicted in the $U$-$J$ phase diagrams in Fig. 3 and Fig. 7 appears to be perpendicular to the two transition lines bordering the metallic phase, i.e., perpendicular to lines with $U-1.51J=\text{constant}$. In Fig. 8, we therefore plot all available data points for the spectral weight $A(0)$ within the HS region of the phase diagram, both for full and reduced band width, as a function of the “effective” interaction parameter $U-1.51J$ (which corresponds to the difference in the local interaction energies between $d^{4}+d^{4}$ and $d^{3}+d^{5}$).

Indeed, all data points plotted in this way essentially condense to a single line. In agreement with the simple atomic-limit prediction (indicated in Fig. 8 as vertical line or shaded regions), the transition to the Mott-insulating phase is shifted towards smaller interaction parameters if the band width is reduced, while the transition to the CD phase is barely affected by this change. However, the occupation difference $\delta n$ in the CD region (not shown) is slightly increased for the case with reduced band width. Furthermore, one can also see that spectral weight disappears only very gradually towards the CDI state, while it disappears rather abruptly at the Mott transition.

We note that, even though the reduction of the $e_{g}$ band width does not significantly affect the CDI phase boundary in our calculations, the band width could nevertheless affect the coupling to the structural distortion and can thus be an important control parameter of the system. In Ref. 20, it was shown that in the two-orbital model the band width controls the relevant electronic susceptibility and is therefore crucial for stabilizing the CD through coupling to the structural distortion.

After establishing and analyzing the phase diagram of the full five-orbital model, we now compare this model to the simplified effective two-orbital description introduced in Sec. II.2. Since by construction the effective two-orbital model assumes a HS configuration on each site, we only compare the HS part of the phase diagram. Our objective is to test how qualitatively and quantitatively accurate the computationally less demanding simplified model is. To ensure a meaningful comparison with the full five-orbital model, the interaction parameters of the two-orbital model are always chosen to be consistent with the corresponding parameters of the full model, as described in Sec. II.2. Furthermore, we only present the case with $\Delta_{e_{g}}=0$, i.e., without an intrinsic $e_{g}$ site splitting between the two sublattices. Our test calculations (not shown) indicate that a site splitting does not alter the quality of the comparison between the full five-orbital model and the effective two-orbital case.

Fig. 9 depicts the occupation difference $\delta n$ and spectral weight $A(0)$ as function of the effective interaction $U-1.51J$, both for the effective two-orbital model and the full five-orbital model. Similar to the case presented in Fig. 8, all points for the same data set, obtained for different combinations of $U$ and $J$, essentially reduce to a single line, demonstrating that $U-1.51J$ is indeed the relevant interaction parameter characterizing the behavior of the model in the HS regime.

Both models give very similar results, with only two small quantitative differences. First, the transition to the CDI phase occurs for slightly smaller effective interaction in the simplified two-orbital model. Due to the large increase of $\delta n$ at the transition, this can lead to quantitative differences in the charge disproportionation of up to 0.5 electrons, but only in a very small region close to the phase transition. Second, the spectral weight at zero energy in the metallic phase is slightly reduced in the effective two-orbital model.

Thus, it appears that the basic observables, such as occupations and the resulting phase boundaries, as well as the gap opening at the Mott transition, obtained from the simplified two-orbital model agree well with the corresponding results for the full five-orbital model. However, this is not true for more complex spectral properties, as can be seen from the $k$-resolved spectral functions shown in Fig. 10 (obtained through analytic continuation of the self-energy).

For the full five-orbital model, the occupied and unoccupied $t_{2g}$ bands appear as rather sharp features. For the case with $U=2.0$ eV and $J=1.0$ eV these bands are situated at around $-$3 eV and 1 eV. The $e_{g}$ spectral weight is much more incoherent, with significant lifetime broadening due to many-body correlations. In the metallic state obtained for $J=1.0$ eV, some dispersion is still visible between $\omega=$-1\text{\,}\mathrm{e}\mathrm{V}$$ and $1\text{\,}\mathrm{e}\mathrm{V}$. At higher energies, a broad incoherent spectral background is present at around 2 eV.

In contrast, the $k$-resolved spectral functions obtained from the effective two-orbital model exhibit much weaker many-body effects. In the metallic regime for $J=1.0$ eV, a sharp band structure is visible, in particular around $\omega=0$. The band dispersion is essentially identical to the one following from the renormalized and spin-mixed hopping amplitudes, also shown as blue lines for comparison in Fig. 10(b). The Coulomb interaction between the $e_{g}$ electrons only leads to a weak lifetime broadening of the individual bands and mostly increases the splitting between local majority and minority spin channels. Close to the CDI transition at $J=2$ eV, the local minority spin bands are shifted to even higher energies and a gap starts opening at zero frequency, consistent with the full five-orbital model. However, in the two-orbital model, the bands around the gap seem to remain rather well defined, while the band edges at approximately $\pm 1$ eV become rather fuzzy, which does not happen in the five-orbital case.

The fact that the effective two-orbital model exhibits much weaker many-body correlations than the full five-orbital model is to be expected, since the interaction with the $t_{2g}$ core spins merely enters via the renormalized hopping amplitudes and in the form of “classical” Zeeman fields. However, the local self-energies also seem to exhibit some non-trivial differences, leading to pronounced lifetime broadening in different energy regions. Nevertheless, these differences do not seem to critically affect integrated quantities such as occupations or the spectral weight around zero energy (see Fig. 9).