Charge self-consistent density functional theory plus ghost rotationally-invariant slave-boson theory for correlated materials

The DFT+gRISB functional is encoded in a Lagrange function Lanatà et al. (2015); Kotliar et al. (2006) represented as follows:

where $N$ is the total number of electrons in the system, determined by the charge-neutrality condition, $\mu$ is the chemical potential, $\rho(\mathbf{r})$ is the electron density, $\mathcal{J}(\mathbf{r})$ is the corresponding constraining field, $E_{\text{Hxc}}$ is the Hartree exchange-correlation functional, $E_{\text{ion-ion}}$ is the ion-ion energy, $E_{\text{ion}}$ is the ionic potential, $E_{dc}^{i}\left[N_{i}^{0}\right]=\frac{U_{i}}{2}N_{i}^{0}\left(N_{i}^{0}-1\right)-\frac{J_{i}}{2}N_{i}^{0}\left(\frac{N_{i}^{0}}{2}-1\right)$ is the double-counting energy functional associated with the $i$-th impurity Czyżyk and Sawatzky (1994); Anisimov et al. (1997b), $N^{0}_{i}$ is the corresponding occupancy, $V_{i}^{0}$ is the corresponding potential, and $U_{i}$ and $J_{i}$ is the corresponding Coulomb interaction and Hund’s coupling interaction, respectively.

The term $\mathcal{L}_{\text{gRISB}}$ is the gRISB Lagrange function associated with following many-body Kohn-Sham-Hubbard “reference system,” expressed in second quantization as follows:

where $x=(\mathbf{r},\sigma)$, $\sigma$ is the spin variable, $\mathbf{r}$ is the position variable, $\int dx$ indicates both the sum over $\sigma$ and the integral over $\mathbf{r}$, $\nabla^{2}$ is the Laplacian, $\hat{\Psi}(x)$ is the Fermionic field operator, $\hat{P}$ is the projector over a generic computational basis span and:

are the annihilation operators associated with the corresponding correlated degrees of freedom, where $\mathbf{R}$ is the unit cell label, $\mathbf{k}$ is the momentum, and $\mathcal{N}$ is the total number of unit cells in the system. The correlated orbital function is denoted by $\phi_{\mathbf{R}i\alpha}(x)$, where $\alpha=1,..,\nu_{i}$ encodes both the orbital degrees of freedom and the spin.

The Lagrangian $\mathcal{L}_{\text{gRISB}}$ can be formally expressed as follows:

where $H^{\text{qp}}$ is the single-particle matrix representation of the so-called quasiparticle Hamiltonian:

for a given computational basis projector $\hat{P}$ and correlated orbital wavefunction $\phi_{\mathbf{R}i\alpha}(x)$, and the local correlated part of the $H^{\text{qp}}$ has the form:

with

The matrix elements of $R_{i}$ and $\lambda_{i}$ are the so-called renormalization coefficients of the quasiparticle Hamiltonian, $\Delta_{i}$ is the quasiparticle single-particle density matrix. The $\hat{{\Psi}}_{\text{u}}(x)$ and $\hat{\Psi}_{\text{c}}$ is the uncorrelated and correlated part of the field operator, respectively, defined as follows:

where $\hat{I}$ is the identity operator, $a=1,..,B\nu_{i}$, with $B$ controlling the accuracy of the gRISB method, and $f_{\mathbf{k}ia}$ are the so-called quasi-particle annihilation operators. We have also introduced the so-called embedding Hamiltonian of the $i$-th impurity:

The matrix elements of $D_{i}$ and $\lambda_{i}^{c}$ is the hybridization and bath coupling constants, respectively, $|\Phi_{i}\rangle$ is the ground state of $\hat{H}^{\text{emb}}_{i}$, $E^{c}_{i}$ is a Lagrange multiplier enforcing the normalization of $|\Phi_{i}\rangle$, and $c_{i\alpha}$, $b_{ia}$ are the impurity and bath Fermionic annihilation operators, respectively. The number of spin orbitals in the bath is $N_{b,i}=B\nu_{i}$.

The charge neutrality is enforced by the chemical potential $\mu$ at quasiparticle occupancy $N+\sum_{i}m_{i}$, where $m_{i}=(B\nu_{i}-\nu_{i})/2$. The reason for this additional term $m_{i}$ is to enforce the physical occupancy to be at the total physical valence number $N$, where the quasiparticle occupancy and the physical occupancy differs by a number $\sum_{i}m_{i}$ Frank et al. (2021); Lanatà (2022). The $\lambda_{i}$ and $\lambda_{i}^{c}$ can also be viewed as the Lagrange multiplier enforcing the gRISB constraints, and $D_{i}$ is a Lagrange multiplier enforcing the structure of the $\left[R_{i}\right]_{,a\alpha}=\langle\Phi_{i}|c^{\dagger}_{i\alpha}b_{ib}|\Phi_{i}\rangle\big{[}\Delta_{i}(1-\Delta_{i})\big{]}_{ba}^{-\frac{1}{2}}$ matrix Lechermann et al. (2007).

The stationary condition of the DFT+gRISB functional leads to the following saddle-point equations:

where $\langle...\rangle_{0}$ denotes the average over the ground state of $H^{\text{qp}}$, and we introduced the following parameterization of the matrices:

where $h_{i,s}$ is an orthonormal basis of the Hermitian matrices Lanatà et al. (2017a). Equation 12 gives rise to the Kohn-Sham potential, and Eq. 14 is the DFT+gRISB charge density, where the local quasiparticle contribution to the density is subtracted and replaced with the contribution from the local physical density matrix Deng et al. (2009). The chemical potential is determined from Eq. 16. The other equations are the standard gRISB equations for model Hamiltonian Lanatà et al. (2017b); Frank et al. (2021); Lanatà (2022); Mejuto-Zaera and Fabrizio (2023); Lee et al. (2023b). The detailed algorithm for solving these equations will be discussed in the next subsection.

With the converged $R_{i}$ and $\lambda_{i}$, one can compute the Green’s function as follows:

where $|0\rangle$ is the vacuum. Equation 24 holds because the $\hat{H}^{\text{qp}}$ is a single-particle Hamiltonian. The spectral function is calculated from $A_{i}(\mathbf{k},\omega)=-\text{Im}G_{i}(\mathbf{k},\omega)/\pi$, and we use a broadening factor of $\eta=0.05$ eV.

The self-energy can be determined from the Dyson equation:

Note that the self-energy is momentum-independent in gRISB, i.e., $\Sigma_{i}(\mathbf{k},\omega)=\Sigma_{i}(\omega)$. The quasiparticle renormalization weight is determined from:

The expectation value of a generic local operator $\hat{O}$ is computed from the embedding wavefunction:

Our implementation closely follows the previous works Lanatà et al. (2015); Lanatà et al. (2017a). We utilize Wien2k for the DFT part of the calculation Blaha et al. (2020). The projector to the correlated orbitals is constructed from the atomic orbital modified from the density functional theory plus embedded dynamical mean-field theory (DFT+eDMFT) code Haule et al. (2010); Yao et al. (2020). The temperature broadening method is utilized for the Brillouin zone integration with a broadening factor of 0.02 eV. The local density approximation functional (LDA) is utilized in our calculation. We use 5000 k-points and 2000 k-points for the NiO and SrVO₃, respectively, and the RK${}_{\text{max}}$ is set to 7. The energy window for constructing the low-energy Hubbard model is [-10 eV — 10e V]. The fully localized limit (FFL) is used as our double-counting scheme Czyżyk and Sawatzky (1994); Anisimov et al. (1997b); Haule (2015), where the nominal valence occupancy is set to 8 and 2 for NiO and SrVO₃, respectively. We use the DMRG approach implemented in the block2 code to solve the ground-state wavefunction of $H^{\text{emb}}$ Zhai et al. (2023). For the DFT+DMFT calculations, we utilize the DFT+eDMFT code with the same parameter setting as in the DFT+gRISB calculations, which provides a consistent benchmark between the DFT+DMFT and DFT+gRISB methods. The continuous-time Quantum Monte Carlo solver is utilized with $10^{9}$ Monte Carlo steps distributed over 200 CPUs, and the temperature is set to $100$ K. For both methods, we treat all five d-orbitals as correlated shells, and the interaction is of full Slater-Condon type Anisimov et al. (1993); Anisimov et al. (1997b).

The DFT+gRISB self-consistent equations are implemented as follows: (1) converge the DFT calculations to obtain the Kohn-Sham eigenvalue and eigenvectors, (2) construct the projector from the Kohn-Sham eigenvector and the local atomic orbitals, (3) solve the gRISB saddle-point equations Eqs. 13-20 with the Kohn-Sham eigenvalues and the projector, (4) use the gRISB saddle-point solution to compute the new charge density from Eq. 14, (4) feedback the new charge density to DFT to update the new exchange-correlated potential (Eq. 12) and go to step (1) until the charge density and total energy is converged. In our calculations, we set the total energy convergence criteria to $10^{-5}$ eV and the charge convergence criteria to $10^{-3}$.

In this section, we apply DFT+gRISB to SrVO₃ and investigate its total energy and electronic structures. This material has been studied extensively in the past decades and serves as an ideal material for benchmarking new approaches Fujimori et al. (1992); Inoue et al. (1995); Yoshimatsu et al. (2010); Liebsch (2003); Pavarini et al. (2004); Sekiyama et al. (2004); Yoshida et al. (2005); Nekrasov et al. (2006); Tomczak et al. (2012); Taranto et al. (2013); Tomczak et al. (2014); Haule et al. (2014); Haule and Birol (2015); Haule (2015); Zhong et al. (2015); Paul and Birol (2019); James et al. (2021); Pickem et al. (2021). It has a cubic perovskite structure, and the main active orbitals around the Fermi level are in the V-$d_{t_{2g}}$ shell. The correlation effect is essential in SrVO₃, leading to significant renormalization of the bandwidth near the Fermi level.

We first discuss the total energy of the DFT+gRISB. Figure 1 summarized the total energy of LDA, LDA+RISB, LDA+gRISB, and LDA+DMFT as a function of the unit cell volume. First, we reproduce the known fact that LDA underestimates the equilibrium volume at $54$ Å, while the experiment observed value is $56.61\ \text{\AA}$. The LDA+RISB improves the equilibrium volume to $56\ \text{\AA}$ towards the experimental value, but the total energy is not consistent with LDA+DMFT at the quantitative level. On the other hand, LDA+gRISB with 15 bath orbitals significantly improves the total energy, in quantitative agreement with the LDA+DMFT results. The equilibrium volume of LDA+gRISB and LDA+DMFT is $56\ \text{\AA}$ and $55.8\ \text{\AA}$, respectively.

We now discuss the electronic structure of SrVO₃. Figure 2 are the momentum-resolved spectral function along the high-symmetry points and the orbital resolved density of state calculated from the LDA, LDA+RISB, LDA+gRISB, and LDA+DMFT approaches at the experimental equilibrium volume. The main characters around the Fermi level are the Vanadium’s $t_{2g}$ orbitals. For the LDA calculation, the bandwidth of the $t_{2g}$ bands is around 2 eV. When including the electronic correlation effects at the LDA+RISB level, we observe slight renormalization of the bandwidth by a factor around $0.8$, which is inconsistent with the LDA+DMFT, where the renormalization factor is around $0.5$. Moreover, the LDA+RISB electronic structure is almost identical to LDA, implying a weakly correlated metal that is inconsistent with LDA+DMFT. This inconsistency can be remedied by utilizing LDA+gRISB with 15 bath orbitals shown in Fig. 2(c). In the LDA+gRISB results, the $t_{2g}$ bands are renormalized by a factor of $0.5$ in agreement with LDA+DMFT. Moreover, the electronic structure closely resembles LDA+DMFT, except for the upper Hubbard band of the $t_{2g}$ orbitals is approximated by coherent bands in LDA+gRISB, which is an established feature of gRISB Lanatà et al. (2017b); Lee et al. (2023a, b).

The total density of states calculated from LDA+RISB, LDA+gRISB, and LDA+DMFT at the equilibrium volume is shown in Fig. 3, which are compared with the photoemission experiment Sawatzky and Allen (1984). The LDA+gRISB and LDA+DMFT density of states around the Fermi level are in good agreement with each other, while the upper Hubbard band in LDA+gRISB shows a coherent feature. All three methods capture the gap between -2 eV and 0.5 eV and the low energy peaks, mainly attributed to the uncorrelated O and Sr atoms. Finally, we show the orbital-resolved quasiparticle weight $Z$ and occupancy in Table 1. The LDA+gRISB and LDA+DMFT values are in quantitative agreement. On the other hand, the LDA+RISB captures reliable occupancy, but the quasiparticle weight is overestimated.

We now apply DFT+gRISB to NiO, which is a prototypical charge-transfer insulator where the Ni-$d_{e_{g}}$ orbitals hybridize with O-$p$ to form the so-called Zhang-Rice band between the lower and the upper Hubbard band(Zaanen et al., 1985; Zhang and Rice, 1988; Kuneš et al., 2007a). The DFT+DMFT approach has been applied to NiO and reveals further insight into electronic structure of this charge-transfer insulator as well as its properties with external pressure, doping, and the surface effects (Kuneš et al., 2007a, b; Zhang et al., 2019; Leonov et al., 2016; Mandal et al., 2019a, b; Leonov et al., 2020; Karolak et al., 2010; Zhu et al., 2020). In this work, we focus on the paramagnetic phase of NiO and compare its total energy and electronic structure with the DFT+DMFT results.

The total energy as a function of the unit cell volume is shown in Fig. 4 for LDA, LDA+RISB, LDA+DMFT, and LDA+DMFT. The LDA significantly underestimates the unit cell volume around 16.6 Å and fails to capture the experimentally observed insulating behavior, which is a well-known feature. The LDA+RISB also fails to produce a Mott insulating solution with realistic Coulomb parameters, $U=10$ eV and $J=1$ eV (Anisimov et al., 1991; Kuneš et al., 2007a), utilized in this work. Therefore, its total energy has a large 1 eV discrepancy compared to LDA+DMFT. On the other hand, LDA+gRISB captures the charge transfer insulator behavior, significantly improving the total energy to a quantitative agreement with the LDA+DMFT values, with a difference of around 10 meV. This difference can be attributed to the finite temperature effect $T=100$ K utilized in the LDA+DMFT calculations.

The DFT bandstructure and density of states are shown in Fig. 5(a). Without breaking the spin symmetry, DFT predicts a metallic solution, which is a known problem in DFT for transition metal oxides (Mattheiss, 1972). The bands around the Fermi level contain the O-$p$ and Ni-$d_{e_{g}}$ orbital. The Ni-$t_{2g}$orbitals are located below the Fermi level and are almost completely filled.

Next, we show the DFT+RISB momentum-resolved spectral function and density of states in Fig. 5(b). We utilize the constrained LDA Coulomb parameters $U=10$ eV and $J=1$ eV in our simulations (Anisimov et al., 1991; Kuneš et al., 2007a). In DFT+RISB, the correlation effects only slightly renormalized the bands around the Fermi level with $Z_{d_{eg}}=0.6$ and $Z_{d_{t2g}}=0.7$ and the system is far from the metal-insulator transition.

The DFT+gRISB momentum-resolved spectral function and density of states are shown in Fig. 5(c). Here, we use 17 bath orbitals in our DFT+gRISB calculations. The density of states shown in Fig. 5(c) demonstrate that DFT+gRISB accurately captures the charge-transfer insulating behavior in the density of state, where the Hubbard bands are opened in the Ni-$d_{e_{g}}$ orbitals and the Zhang-Rice band is observed around $-1$ eV with strong $O$-$p$ and $Ni$-$d$ hybridization, in good agreement with the DFT+DMFT results shown in Fig. 5(d) and the previous studies (Kuneš et al., 2007a, b; Zhang et al., 2019; Leonov et al., 2016, 2020; Karolak et al., 2010; Zhu et al., 2020). Moreover, the DFT+gRISB momentum-resolved spectral functions $A(\mathbf{k},\omega)$ captures reliably the dispersive excitations compared to DFT+DMFT, except the incoherent broadening features, which cannot be described within the gRISB framework.

In Fig. 6, we compare the LDA+gRISB total density of states with LDA+DMFT and the photoemission bremsstrahlung-isochromat-spectroscopy (XPS/BIS) (Sawatzky and Allen, 1984). Our LDA+gRISB density of states captures the main features in the XPS/BIS spectrum, where the band gap is about $4$ eV, and the heights of the peaks on the band edges are reasonably captured. On the other hand, LDA+DMFT has band gap around 2 eV, smaller than the experimental band gap.

Finally, the orbital-resolved occupancy of different approaches is shown in Table 2. The LDA+gRISB’s occupancy is in good agreement with LDA+DMFT and improves the LDA values. On the other hand, although LDA+RISB fails to capture the charge-transfer insulating solution, its occupancy is identical to the LDA+gRISB values and close to the LDA+DMFT values.