Emergent Bloch Excitations in Mott Matter

TL;DR

This paper presents a unified theoretical framework for understanding excitations in Mott insulators, revealing emergent Bloch-like features and a topological transition at the Mott point, with implications for modeling strongly correlated materials.

Contribution

It introduces a non-perturbative, unified picture of excitations in Mott systems, connecting quasiparticles and Hubbard bands within an emergent Fermi liquid framework.

Findings

Incoherent excitations exhibit well-defined Bloch character.

The Mott transition is a topological change in dispersing bands.

The variational approach yields quantitatively accurate results.

Abstract

We develop a unified theoretical picture for excitations in Mott systems, portraying both the heavy quasiparticle excitations and the Hubbard bands as features of an emergent Fermi liquid state formed in an extended Hilbert space, which is non-perturbatively connected to the physical system. This observation sheds light on the fact that even the incoherent excitations in strongly correlated matter often display a well defined Bloch character, with pronounced momentum dispersion. Furthermore, it indicates that the Mott point can be viewed as a topological transition, where the number of distinct dispersing bands displays a sudden change at the critical point. Our results, obtained from an appropriate variational principle, display also remarkable quantitative accuracy. This opens an exciting avenue for fast realistic modeling of strongly correlated materials.

Figures

Equations

$$\hat{H}$$

$$⟨{\Psi }_0|{\hat{\mathcal{P}}}_R^{†}{\hat{\mathcal{P}}}_R^{}|{\Psi }_0⟩$$

$$⟨{\Psi }_0|{\hat{\mathcal{P}}}_R^{†}{\hat{\mathcal{P}}}_R^{}{c}_{R\alpha \sigma }^{†}{c}_{R\beta \sigma }^{}|{\Psi }_0⟩$$

$${\hat{H}}_{\text{qp}}$$

$$G(\epsilon_{k},\omega)\!=\!\left[\tilde{\mathcal{R}}^{\dagger}\frac{1}{\omega\!-\!\big{(}\tilde{\mathcal{R}}\tilde{\epsilon}_{k}\tilde{\mathcal{R}}^{\dagger}\!+\!\tilde{\lambda}\big{)}}\tilde{\mathcal{R}}\right]_{11}\!\!\!\!=\!\big{[}\omega\!-\!\epsilon_{k}\!-\!\Sigma(\omega)\big{]}^{-1}\!\!,$$

$${\tilde{\lambda }}_{ij}$$

$${\tilde{\mathcal{R}}}_{ij}$$

$$\Sigma (\omega )=\frac{\omega }{1+\frac{1}{z-\frac{{\omega }^2-l^2+2h{\omega }^2}{{\omega }^2-l^2}}}=-\frac{1-z}{z}\omega +o({\omega }^2).$$

$${\tilde{\lambda }}_{ij}$$

$${\tilde{\mathcal{R}}}_{ij}$$

$$\Sigma (\omega )=-\frac{1-h}{h}\omega +\frac{l^2}{h}\frac{1}{\omega }.$$

$$\hat{H}=\sum _{R{R}^{\prime }}\sum _{\alpha ,\beta =1}^{\tilde{\nu }}{t}_{R{R}^{\prime }}^{\alpha \beta }{c}_{R\alpha }^{†}{c}_{R^{\prime }\beta }^{}+\sum _R{\hat{H}}_R^{\text{loc}}=\sum _k\sum _{\alpha ,\beta =1}^{\tilde{\nu }}{ϵ}_k^{\alpha \beta }{c}_{k\alpha }^{†}{c}_{k\beta }^{}+\sum _R{\hat{H}}_R^{\text{loc}},$$

$$|{\Psi }_G⟩={\hat{\mathcal{P}}}_G|{\Psi }_0⟩,$$

$${\hat{\mathcal{P}}}_R=\sum _{A,n=0}^{2^{\tilde{\nu }}-1}{\tilde{\Lambda }}_{An}|A,R⟩⟨n,R|,$$

$$|A,R⟩$$

$$⟨n,R|$$

$$⟨{\Psi }_0|f_{Ra}^{†}{f}_{Rb}^{}|{\Psi }_0⟩\propto {\delta }_{ab}\forall a,b\in \{1,...,\tilde{\nu }\}.$$

$$⟨{\Psi }_0|{\hat{\mathcal{P}}}_R^{†}{\hat{\mathcal{P}}}_R^{}|{\Psi }_0⟩$$

$$⟨{\Psi }_0|{\hat{\mathcal{P}}}_R^{†}{\hat{\mathcal{P}}}_R^{}{f}_{R\alpha }^{†}{f}_{R\beta }^{}|{\Psi }_0⟩$$

$$\mathcal{L}[{\Psi }_0,E;\Phi ,E^c;\mathcal{R},{\mathcal{R}}^{†},\lambda ;\mathcal{D},{\mathcal{D}}^{†},{\lambda }^c;\Delta ]=\frac{1}{\mathcal{N}}⟨{\Psi }_0|{\hat{H}}_{\text{qp}}[\mathcal{R},{\mathcal{R}}^{†};\lambda ]|{\Psi }_0⟩+E\left(1-⟨{\Psi }_0|{\Psi }_0⟩\right)$$

$$+\left[⟨\Phi |{\hat{H}}_{\text{emb}}[\mathcal{D},{\mathcal{D}}^{†};{\lambda }^c]|\Phi ⟩+E^c\left(1-⟨\Phi |\Phi ⟩\right)\right]$$

$$-\left[\sum _{a,b=1}^{\tilde{\nu }}\left({\lambda }_{ab}+{\lambda }_{ab}^c\right){\Delta }_{ab}+\sum _{c,a,\alpha =1}^{\tilde{\nu }}\left({\mathcal{D}}_{a\alpha }{\mathcal{R}}_{c\alpha }{\left[\Delta (1-\Delta )\right]}_{ca}^{\frac{1}{2}}+\text{c.c.}\right)\right],$$

$${\hat{H}}_{\text{qp}}$$

$${\hat{H}}_{\text{emb}}$$

$$|\Phi ⟩\equiv \sum _{A,n=0}^{2^{\tilde{\nu }}-1}{e}^{i\frac{\pi }{2}{N}_n(N_n-1)}{\tilde{\varphi }}_{An}[{{\hat{c}}^{†}}_1{]}^{q_1(A)}...[{{\hat{c}}^{†}}_{\tilde{\nu }}{]}^{q_{\tilde{\nu }}(A)}[{{\hat{f}}^{}}_1{]}^{q_1(n)}...[{{\hat{f}}^{}}_{\tilde{\nu }}{]}^{q_{\tilde{\nu }}(n)}|\bar{0}⟩,$$

$$\tilde{\varphi }=\tilde{\Lambda }\sqrt{P^0},$$

$$P_{nm}^0=⟨{\Psi }_0|[f_{R\tilde{\nu }}^{}{]}^{q_{\tilde{\nu }}(n)}...[f_{R1}^{}{]}^{q_1(n)}[f_{R1}^{†}{]}^{q_1(m)}...[f_{R\tilde{\nu }}^{†}{]}^{q_{\tilde{\nu }}(m)}|{\Psi }_0⟩,$$

$$\Delta$$

$${\lambda }^c$$

$$\lambda$$

$$\mathcal{R}$$

$$\displaystyle\frac{1}{\mathcal{N}}\sum_{k}\big{[}f\!\left(\mathcal{R}\epsilon_{k}\mathcal{R}^{\dagger}+\lambda\right)\big{]}_{ba}=\Delta_{ab}$$

$$\displaystyle\frac{1}{\mathcal{N}}\,\sum_{k}\big{[}\epsilon_{k}\mathcal{R}^{\dagger}\,f\!\left(\mathcal{R}\epsilon_{k}\mathcal{R}^{\dagger}+\lambda\right)\big{]}_{\alpha a}=\sum_{c}\mathcal{D}_{c\alpha}\left[\Delta\left(1-\Delta\right)\right]^{\frac{1}{2}}_{ac}$$

$$\sum _{cb\alpha }\frac{\partial }{\partial d_s^p}{\left[\Delta \left(1-\Delta \right)\right]}_{cb}^{\frac{1}{2}}{\mathcal{D}}_{b\alpha }{\mathcal{R}}_{c\alpha }+\text{c.c.}+{\left[l+l^c\right]}_s=0$$

$${\hat{H}}^{\text{emb}}[\mathcal{D},{\lambda }^c]|\Phi ⟩=E^c|\Phi ⟩$$

$$\displaystyle\big{[}\mathcal{F}^{(1)}\big{]}_{\alpha a}\equiv\langle\Phi|\,\hat{c}^{\dagger}_{\alpha}\hat{f}^{\phantom{\dagger}}_{a}\,|\Phi\rangle-\sum_{c}\left[\Delta\left(1-\Delta\right)\right]^{\frac{1}{2}}_{ca}[\mathcal{R}]_{c\alpha}=0$$

$$\displaystyle\big{[}\mathcal{F}^{(2)}\big{]}_{ab}\equiv\langle\Phi|\,\hat{f}^{\phantom{\dagger}}_{b}\hat{f}^{\dagger}_{a}\,|\Phi\rangle-\left[\Delta\right]_{ab}=0\,,$$

Full text

Introduction:— The physical nature of the excited states in strongly interacting quantum systems has long been a subject of much controversy and debate. Deeper understanding was achieved by Landau, more than half a century ago Landau (1957), who realized that in systems of fermions the Pauli principle provides a spectacular simplification. He showed that many properties of Fermi systems can be understood in terms of weakly interacting quasiparticles (QP), allowing a precise and detailed description of strongly correlated matter. Modern experiments provide for even more direct evidence of such QP excitations, for example from using angle-resolved photoemission spectroscopy (ARPES) Damascelli (2004) or scanning-tunneling microscopy (STM) methods Binnig et al. (1982).

The Fermi liquid paradigm, however, describes only the low energy excitations. At higher energies, the physical properties are often dominated by incoherent processes, which do not conform to the Landau picture. The task to provide a simple and robust theoretical description of such incoherent excitations has therefore emerged as a central challenge of contemporary physics. An intriguing apparent paradox is most evident around the Mott point. Here, ARPES and STM experiments provide often clear evidence of additional well-defined high energy excitations (Hubbard bands) which, while being fairly incoherent, still display relatively well defined Bloch character with pronounced momentum dispersion, see, e.g., Ref. Evtushinsky et al. (2017). As a matter of fact, it is often difficult to experimentally even distinguish the Hubbard bands found in Mott insulators from ordinary Bloch bands found at high energy in conventional band insulators. While such behavior can be already numerically reproduced by some modern many-body approximations Georges et al. (1996); Anisimov and Izyumov (2010), a simple conceptual picture for the apparent Bloch character of such high energy charge excitations is not still available. In particular, variational methods such as the Gutzwiller Approximation (GA) Gutzwiller (1965) — which are often able to reproduce the numerical results in a much simpler semi-analytical fashion — generally capture only the low-lying QP features on the metallic side, but cannot provide a description of charge excitations around the Mott point and in the insulating regime.

The goal of this Letter is to write an appropriate variational wave function able to capture both the (low energy) QP bands and the (high energy) Hubbard bands, within the same theoretical framework. A particularly interesting fact emerging from our theory is that many important features of both types of excitations are encoded in the bare density of states (DOS) of the uncorrelated system and a few renormalization parameters — in a similar fashion as for the QP excitations in Landau theory of Fermi liquids. This is accomplished, similarly as in many other theories for many-body systems, see, e.g., Refs. Affleck et al. (1987); Vaestrate et al. (2008); Schollwock (2011), by enlarging the Hilbert space by introducing auxiliary ”ghost” degrees of freedom. In particular, this construction sheds light on the physical origin of the ”hidden” Bloch character of the Hubbard bands mentioned above. Our calculations of the single-band Hubbard model, which are benchmarked against the Dynamical Mean Field Theory (DMFT) Georges et al. (1996); Anisimov and Izyumov (2010) solution, show that the new wave function quantitatively captures not only the dispersion of the QP but also the Hubbard bands. Furthermore, our theory enables us to describe the Mott transition and the coexistence region between the metallic and the Mott-insulator phases.

Ghost GA theory:— For simplicity, our theory will be formulated here for the single-band Hubbard model

[TABLE]

at half-filling. The generalization to arbitrary multi-orbital Hubbard Hamiltonians is straightforward sup .

In order to construct the Ghost-GA theory we are going to embed the physical Hamiltonian of the system [Eq. (1)] within an extended Hilbert space obtained by introducing auxiliary Fermionic ”ghost” degrees of freedom not coupled with the physical orbitals, see Fig. 1.

Let us represent $\hat{H}$ within the extended Hilbert space mentioned above as follows:

[TABLE]

where $\tilde{t}_{RR^{\prime}}^{11}=t_{RR^{\prime}}$ are the physical hopping parameters, $\tilde{\epsilon}_{k}^{11}=\epsilon_{k}$ are the eigenvalues of the first term of $\hat{H}$, $\tilde{t}_{RR^{\prime}}^{\alpha\beta}=\tilde{\epsilon}_{k}^{\alpha\beta}=0$ $\forall\,(\alpha,\beta)\neq(1,1)$ and $\sigma$ is the spin.

Our theory consists in applying the ordinary multi-orbital GA theory Lanatà et al. (2008); Bünemann et al. (1998); Lanatà et al. (2015a); Attaccalite and Fabrizio (2003); Deng et al. (2009) to Eq (2). In other words, the expectation value of $\hat{H}$ is optimized variationally with respect to a Gutzwiller wave function represented as $|\Psi_{G}\rangle=\hat{\mathcal{P}}_{G}|\Psi_{0}\rangle$, where $|\Psi_{0}\rangle$ is the most general Slater determinant, $\hat{\mathcal{P}}_{G}=\prod_{R}\hat{\mathcal{P}}_{R}$, and $\hat{\mathcal{P}}_{R}$ acts over all of the local degrees of freedom labeled by $R$ — including the ghost orbitals $\alpha>1$. The variational wave function is restricted by the following conditions:

[TABLE]

which are commonly called ”Gutzwiller constraints”. Furthermore, the so called ”Gutzwiller Approximation” Gutzwiller (1965) — which is exact in the limit of infinite dimensions (where DMFT is exact) — is employed. The minimization of the variational energy will be performed by employing the algorithms derived in Ref. Lanatà et al. (2017).

The basis of our theory is that extending the Hilbert space by introducing the ghost orbitals does not affect the physical Hubbard Hamiltonian $\hat{H}$, as all of its terms involving ghost orbitals are multiplied by [math], see Eq. (2). The advantage of enlarging the Hilbert space arises exclusively from the fact that the corresponding Ghost-GA multi-orbital variational space is substantially more rich with respect to the original GA variational space (where $\hat{\mathcal{P}}_{R}$ acts only on the physical degrees of freedom) sup . In this respect, our scheme presents analogies, e.g., with the solution of the Affleck, Lieb, Kennedy and Tasak (AKLT) model Affleck et al. (1987) and with the theories of Matrix Product States (MPS) and Projected Entangled Pair States (PEPS) Vaestrate et al. (2008); Schollwock (2011), which are also variational constructions involving virtual entanglement and local maps foo (a). Further technical details about the role of the ghost orbitals are provided in the supplemental material sup .

As shown in previous works, see, e.g., Refs. Lanatà et al. (2015a), the variational energy minimum of $\hat{H}$ is realized by a wave function $|\Psi_{G}\rangle=\hat{\mathcal{P}}_{G}|\Psi_{0}\rangle$ where $|\Psi_{0}\rangle$ is the ground state of a quadratic multi band Hamiltonian represented as

[TABLE]

where $f_{ka\sigma}$ are related to $c_{ka\sigma}$ by a proper unitary transformation Lanatà et al. (2008); Lanatà et al. (2015a), the matrices $\tilde{\mathcal{R}}$ and $\tilde{\lambda}$ are determined variationally and $\tilde{\epsilon}^{*}_{kn}$ are the eigenvalues of $\hat{H}_{\text{qp}}$. The states $|\Psi^{p}_{Gkn\sigma}\rangle=\hat{\mathcal{P}}_{G}\psi^{\dagger}_{kn\sigma}|\Psi_{0}\rangle$ and $|\Psi^{h}_{Gkn\sigma}\rangle=\hat{\mathcal{P}}_{G}\psi^{\phantom{\dagger}}_{kn\sigma}|\Psi_{0}\rangle,$ where $\psi^{\dagger}_{ka\sigma}$ are the eigen-operators of $\hat{H}_{\text{qp}}$, represent excited states of $\hat{H}$ Bünemann et al. (2003); Lanatà (2009); sup .

The energy-resolved Green’s function of the physical degrees of freedom ($\alpha=1$) can be evaluated in terms of the excitations mentioned above sup and represented as

[TABLE]

where the subscript ”$11$” indicates that we are interested only in the physical component $\alpha=\beta=1$ of the Green’s function. We point out that, since Eq. (6) involves a matrix inversion Savrasov et al. (2006, 2005), the Ghost-GA approximation $\Sigma(\omega)$ to the physical self-energy is generally a non-linear function sup — while it is linear by construction in the ordinary GA theory. Note also that the poles of $G(\epsilon_{k},\omega)$ coincide with the eigenvalues $\tilde{\epsilon}_{kn}^{*}$ of $\hat{H}_{\text{qp}}$, see Eq. (5).

Application to the single-band Hubbard model:— Below we apply our approach to the Hubbard Hamiltonian [Eq. (2)] at half-filling assuming a semicircular DOS foo (b), which corresponds, e.g., to the Bethe lattice in the limit of infinite connectivity, where DMFT is exact Georges et al. (1996). The half-bandwidth $D$ will be used as the unit of energy. The extended Ghost-GA scheme will be applied following the procedure of Ref. Lanatà et al. (2017), utilizing up to 2 ghost orbitals.

In Fig. 2 is shown the evolution as a function of the Hubbard interaction strength $U$ of the Ghost-GA total energy, the local double occupancy and the QP weight $z$. Our results are shown in comparison with the ordinary GA theory and with DMFT in combination with Numerical Renormalization Group (NRG). In particular, we employed the ”NRG Ljubljana” impurity solver Zitko (2011).

The agreement between Ghost-GA and DMFT is quantitatively remarkable. In particular, the Ghost-GA theory enables us to account for the coexistence region of the Mott and metallic phases, which is not captured by the ordinary GA theory. The values of the boundaries of the coexistence region $U_{c1}\simeq 2$, $U_{c2}\simeq 2.88$ are in good agreement with the DMFT results available in the literature García et al. (2004); Werner and Millis (2007); Bulla (1999); Tong et al. (2001), i.e., $U_{c1}\simeq 2.39$, $U_{c2}\simeq 2.94$. The Ghost-GA value of $U_{c2}$, which is the actual Mott transition point at $T=0$, is particularly accurate. The method also gives a reasonable value for the very small energy scale characterizing the coexistence region, which we can estimate as $T_{c}\simeq E_{\text{ins}}(U_{c1})-E_{\text{met}}(U_{c1})\simeq 0.02$, consistently with both DMFT and experiments Aguiar et al. (2005); Limelette et al. (2003). We point out also that, as shown in the second panel of Fig. 2, the Ghost-GA approach captures the charge fluctuations in the Mott phase, while this is approximated by the simple atomic limit (which has zero double occupancy) within the Brinkman-Rice scenario Brinkman and Rice (1970).

Interestingly, while at least 2 ghost orbitals are necessary to obtain the data illustrated above for the Metallic solution, 1 ghost orbital is sufficient to obtain our results concerning the Mott phase. Increasing further the number of ghost orbitals does not lead to any appreciable difference sup . As we are going to show, this is connected with the fact that the electronic structures of the Mott and the metallic phases are topologically distinct.

Let us now analyze the Ghost-GA single-particle Green’s function $G(\epsilon,\omega)$, see Eq. (6). In Fig. 3 is shown the Ghost-GA energy-resolved spectral function $A(\epsilon,\omega)=-\frac{1}{\pi}\text{Im}{G(\epsilon,\omega)}$ in comparison with DMFT Karski et al. (2008).

Although the broadening of the bands (scattering rate), is not captured by our approximation (as it is not captured by the ordinary GA), the positions and the weights of the poles of the Ghost-GA spectral function encode most of the DMFT features, not only at low energies (QP excitations), but also at high energies (Hubbard bands). In order to analyze how the spectral properties of the system emerge within the Ghost-GA theory, it is particularly convenient to express the QP Hamiltonian [Eq. (5)] in a gauge where $\tilde{\lambda}$ is diagonal foo (c).

In the metallic phase, an explicit Ghost-GA calculation obtained employing 2 ghost orbitals shows that the matrices $\tilde{\mathcal{R}}$ and $\tilde{\lambda}$ are represented as follows:

[TABLE]

where $\delta_{ij}$ is the Kronecker delta, and $l$, $z$ and $h$ are real positive numbers determined numerically as in Ref. Lanatà et al. (2017). The corresponding self-energy, see Eq. (6), is foo (d):

[TABLE]

Thus, the variational parameter $z$ of Eq. (8) represents the QP weight, whose behavior was displayed in the third panel of Fig. 2. Note that the overall spectral weight $\int d\omega\int d\epsilon\,\rho(\epsilon)\,A(\epsilon,\omega)$, where $\rho(\epsilon)$ is the semicircular DOS, is not $z$ as in the ordinary GA theory, but it is $z+h=[\tilde{\mathcal{R}}^{\dagger}\tilde{\mathcal{R}}]_{11}$, which is almost equal to $1$ for all values of $U$ (see the inset of the third panel in Fig. 2). The additional spectral contribution $h$, which is not present in the ordinary GA approximation, enables the Ghost-GA theory to account for the Hubbard bands.

In the Mott phase, an explicit Ghost-GA calculation obtained employing 1 ghost orbital shows that the matrices $\tilde{\mathcal{R}}$ and $\tilde{\lambda}$ are represented as follows:

[TABLE]

where $l$ and $h$ are real positive numbers determined numerically as in Ref. Lanatà et al., 2017. Note that $h=[\tilde{\mathcal{R}}^{\dagger}\tilde{\mathcal{R}}]_{11}\simeq 1$ (see the inset of the third panel in Fig. 2). The corresponding self-energy, see Eq. (6), is foo (d):

[TABLE]

The pole of the self energy at $\omega=0$, which is the source of the Mott gap, is captured by the Ghost-GA theory.

The analysis above clarifies also why, by construction, within the Ghost-GA approximation the self-energy can develop poles, see Eqs. (9) and (12), but can not capture branch-cut singularities on the real axis.

It is important to point out that the Hilbert space extension that we have introduced in this work has been essential in order to capture the effect of the electron correlations on the topology of the excitations — such as the change of the number of bands at the Mott transition (between 3 bands in the metallic phase and 2 bands in the Mott phase). In fact, without extending the Hilbert space, the ordinary GA theory enables only to renormalize and shift the band structure with respect to the uncorrelated limit $U=0$, without affecting its qualitative topological structure. On the other hand, extending the Hilbert space enables us to relax this constraint, as $G(\epsilon,\omega)$, see Eq. (6), is variationally allowed to have any number of distinct poles equal or smaller to the corresponding total (physical and ghost) number of orbitals sup . It is for this reasons that only 1 ghost (2 orbitals) is sufficient to describe the Mott phase of the single-band Hubbard model, while at least 2 ghosts (3 orbitals) are necessary in order to describe its metallic phase — whose spectra includes the QP excitations and the 2 Hubbard bands. A remarkable aspect of this construction is that, within the Ghost-GA theory, the information concerning the spectral function — including the Hubbard bands — is entirely encoded in only 3 parameters ($z,h,l$) in the Metallic phase, and in 2 parameters ($h,l$) in the Mott phase, see Eqs. (7), (8), (10), (11).

Conclusions:— We derived a unified theoretical picture for excitations in Mott systems based on a generalization of the GA, which captures not only the low-energy QP excitations, but also the Hubbard bands. The essential idea consists in extending the Hilbert space of the system by introducing auxiliary ”ghost” orbitals. This construction enables us to express analytically many important features of both types of excitations in terms of the bare DOS of the uncorrelated system and a few renormalization parameters, in a similar fashion as for the QP excitations in Landau theory of Fermi liquids. In particular, this idea provides us with a conceptual picture which assigns naturally a Bloch character to the Hubbard bands even in Mott insulators. In this respect, we note that our theory presents a few suggestive analogies with the interesting (but rather speculative) idea of ”hidden Fermi liquid” previously introduced by P. W. Anderson Anderson (2008) within the context of the BCS wave function (for superconductors) and the Laughlin’s Jastrow wave function (for the Fractional Hall Effect). In fact, they both propose a descriptions of non-Fermi liquid states related to ordinary Fermi liquids residing in unphysical Hilbert spaces, see, e.g., Ref. Jain and Anderson (2009). From the computational perspective, our Ghost-GA theory constitutes a very promising tool for ab-initio calculations in combination with Density Functional Theory (DFT) Hohenberg and Kohn (1964); Kohn and Sham (1965); Deng et al. (2009); Ho et al. (2008); Lanatà et al. (2015a), as it is substantially more accurate with respect to the ordinary GA approximation, without much additional computational cost. In fact, within the numerical scheme described in Refs. Lanatà et al. (2017); Lanatà et al. (2015a), our theory results in solving iteratively a finite impurity model, where the number of bath sites grows linearly with the total number of ghost orbitals sup . Since there exist numerous available techniques enabling to solve efficiently this auxiliary problem, see, e.g., Refs. Wolf et al. (2015); Saberi et al. (2008); Weichselbaum et al. (2009); Lu et al. (2014); Lanatà et al. (2016), this opens an exciting avenue for realistic modeling of many challenging materials, including predictions of ARPES spectra for complex orbitally-selective Mott insulators. Furthermore, since the Ghost-GA theory is based on the multi-orbital GA Lanatà et al. (2017); Lanatà et al. (2015a), it can be straightforwardly generalized to finite temperatures Lanatà et al. (2015b); Sandri et al. (2013); Wang et al. (2010), to non-equilibrium problems Schirò and Fabrizio (2010); Lanatà and Strand (2012), and to calculate linear response functions Fabrizio (2017). For the same reason, the Ghost-GA theory can be also straightforwardly reformulated Bünemann and Gebhard (2007); Lanatà et al. (2008); Lanatà et al. (2015a) in terms of the rotationally invariant slave boson (RISB) theory Lechermann et al. (2007); Li et al. (1989); sup , whose exact operatorial foundation recently derived in Ref. Lanatà et al. (2017) constitutes a starting point to calculate further corrections Dao and Frésard (2017). It would be also interesting to apply the ghost-orbital Hilbert space extension in combination with the Variational Monte Carlo method Ceperley et al. (1977) or the generalization of the GA to finite dimensions of Ref. Münster and Bünemann (2016), which might lead to a more accurate description of strongly correlated electron systems even beyond the DMFT approximation.

1 I.   Technical details about Ghost-GA theory

1.1 A.   Introduction

Let us consider a generic multi-orbital Hubbard Hamiltonian represented as

[TABLE]

where $\tilde{\nu}$ is the total number of Fermionic degrees of freedom (both spin and orbital), and $\hat{H}_{R}^{\text{loc}}$ is a generic local operator (i.e., any operator that can be expressed in terms of operators $\{c^{\dagger}_{R\alpha},c^{\phantom{\dagger}}_{R\beta}\}$ with fixed label $R$).

The multi-orbital Gutzwiller Approximation [Gutzwiller, 1965; Lanatà et al., 2008; Bünemann et al., 1998; Lanatà et al., 2015; Attaccalite and Fabrizio, 2003] (GA) consists in minimizing the expectation value of $\hat{H}$ with respect the most general Gutzwiller wave function:

[TABLE]

where $|\Psi_{0}\rangle$ is the most general Slater determinant, $\hat{\mathcal{P}}_{G}=\prod_{R}\hat{\mathcal{P}}_{R}$, and $\hat{\mathcal{P}}_{R}$ is the most general map acting over all of the local degrees of freedom labeled by $R$. As discussed in Ref. [Lanatà et al., 2015], it is convenient to represent $\hat{\mathcal{P}}_{R}$ in a ”mixed basis” representation as follows:

[TABLE]

where

[TABLE]

the occupation numbers $q_{1}(l),...,q_{\tilde{\nu}}(l)$ are the digits of a generic integer $l\in\{0,...,2^{\tilde{\nu}}-1\}$ in basis $2$, and the ladder operators $f^{\dagger}_{Ra}$ are related to $c^{\dagger}_{R\alpha}$ through a unitary transformation such that

[TABLE]

The variational wave function is restricted by the following conditions:

[TABLE]

which are commonly called ”Gutzwiller constraints”. Furthermore, the ”Gutzwiller Approximation”, which becomes exact in the limit of infinite dimension [Gutzwiller, 1965] — where Dynamical Mean Field Theory (DMFT) is exact [Georges et al., 1996] — is assumed.

As shown in previous works, see, e.g., Refs. [Lanatà et al., 2008, 2016], the accuracy of the Gutzwiller method is improved substantially when the operator $\hat{\mathcal{P}}_{G}$ (Gutzwiller projector) is allowed to act also beyond the space generated by the correlated degrees of freedom. The main technical idea underlying the Ghost-GA theory presented in this work consists essentially in creating new auxiliary (inert) orbitals with the sole purpose of enriching the variational space without changing its formal structure (i.e., a Gutzwiller-projected Slater determinant). Of course, as the ordinary GA, the Ghost-GA theory exploits the so called ”Gutzwiller Approximation” [Gutzwiller, 1965], which becomes exact in the limit of infinite dimension (where DMFT is exact). Thus, in summary, the Ghost-GA theory consists in an improved variational approximation to DMFT realized by enlarging the ordinary GA variational space.

Further details about the role of the ghost orbitals are provided in the subsections below.

1.2 B.   Functional formulation of the Gutzwiller energy minimization

As discussed in the main text, all of our calculations of the single band Hubbard Hamiltonian [see Eq. 1 of the main text] have been performed by applying the ordinary multi-orbital Gutzwiller theory to the same Hamiltonian ”embedded” within the extended Hilbert space obtained by introducing auxiliary ghost orbitals [see Eq. 2 of the main text]. Specifically, we have employed the numerical procedure derived in described in Refs. [Lanatà et al., 2017, 2015], which we summarize below for completeness.

For simplicity, the theory will be formulated here assuming a translationally invariant Gutzwiller solution. As shown in Refs. [Lanatà et al., 2017, 2015], the minimization of the expectation value of a generic multi-orbital Hamiltonian represented as in Eq. (1) with respect to the Gutzwiller wave function can be conveniently formulated in terms of the following Lagrange function:

[TABLE]

where $\mathcal{N}$ is the total number of $k$ points, $E$ and $E^{c}$ are real numbers, $\Delta$, $\lambda^{c}$ and $\lambda$ are Hermitian matrices, $\mathcal{D}$ and $\mathcal{R}$ are generic complex matrices. The auxiliary Hamiltonians $\hat{H}_{\text{qp}}$ and $\hat{H}_{\text{emb}}$, which are called ”quasiparticle Hamiltonian” and ”Embedding Hamiltonian”, respectively, are defined as follows:

[TABLE]

The auxiliary Fermionic Hamiltonian $\hat{H}_{\text{emb}}$, which was introduced in Ref. [Lanatà et al., 2015], enables us to exploit techniques such as those developed in quantum chemistry in order to tackle the energy minimization problem. The state $|\Phi\rangle$, see Eq. (9), is the most general many-body state — with $\tilde{\nu}$ electrons — belonging to the auxiliary embedding Hilbert space mentioned above.

1.3 C.   Connection between embedding state $|\Phi\rangle$ and Gutzwiller variational parameters

Let us consider the following expansion of the most general state $|\Phi\rangle$:

[TABLE]

where $A,n\in\{0,...,2^{\tilde{\nu}}-1\}$, $N_{n}=\sum_{a=1}^{\tilde{\nu}}q_{a}(n)$ and $|\bar{0}\rangle=[{\hat{f}^{\dagger}}_{1}]\,...\;[{\hat{f}^{\dagger}}_{\tilde{\nu}}]\,|0\rangle$.

For later convenience, it is important to point out that, as shown in Ref. [Lanatà et al., 2015], the matrix $\tilde{\phi}$ with entries $\tilde{\phi}_{An}$ is connected with the Gutzwiller variational parameters as follows:

[TABLE]

where $\tilde{\Lambda}$ is the matrix defining $\hat{\mathcal{P}}_{R}$ as in Eq. (3), $P^{0}$ is $2^{\tilde{\nu}}\times 2^{\tilde{\nu}}$ matrix with entries

[TABLE]

and $n,m\in\{0,...,2^{\tilde{\nu}}-1\}$. Note that, because of Eq. (6), $P^{0}$ is diagonal.

1.4 D.   Gutzwiller Lagrange equations in the standard multi-orbital GA

Following Ref. [Lanatà et al., 2015], in order to take into account the fact that $\Delta$, $\lambda^{c}$ and $\lambda$ are Hermitian matrices we introduce the following parametrizations:

[TABLE]

where the set of matrices $h_{s}$ is an orthonormal basis of the space of Hermitian matrices with dimension $\nu$, ${}^{t}h_{s}$ are the corresponding transposed matrices, and $d_{s}$, $l_{s}^{c}$ and $l_{s}$ are real numbers, while $r_{s}$ are complex numbers. The above-mentioned orthonormality is defined with respect to the standard scalar product $(A,B)\equiv\mathrm{Tr}\!\left[A^{\dagger}B\right]$.

It can be readily shown that the saddle-point conditions of $\mathcal{L}$, see Eq. (9), with respect to all of its arguments provides the following system of Lagrange equations:

[TABLE]

where the function $f$ appearing in Eqs. (19) and (20) is the Fermi function at zero temperature ($T=0$).

The solution of the equations above can be obtained numerically employing the following procedure [Lanatà et al., 2015]. (I) Given a set of coefficients $r_{s}$ and $l_{s}$, determine the corresponding matrices $\mathcal{R}$ and $\lambda$ using Eqs. (17) and (18), and calculate $\Delta$ using Eq. (19). (II) Calculate $\mathcal{D}$ by inverting Eq. (20). (III) Calculate the coefficients $l^{c}_{s}$ using Eq. (21) and the corresponding matrix $\lambda^{c}$ using Eq. (16). (IV) Construct the embedding Hamiltonian $\hat{H}^{\text{emb}}$ and compute its ground state $|\Phi\rangle$, see Eq. (22), within the subspace with $\nu$ electrons. (V) Determine the left members of Eqs. (23) and (24). The equations (23) and (24) are satisfied if and only if the coefficients $r_{s}$ and $l_{s}$ proposed at the first of the steps above identify a solution of the GA Lagrange function. The steps above enable us to reduce the solution of the GA Lagrange equations to a root problem, which can be formally represented as

[TABLE]

that can be readily solved numerically, e.g., using the quasi-Newton method.

All of the expectation values with respect to the Gutzwiller wavefunction can be readily expressed in terms of the states $|\Phi\rangle$ and $|\Psi_{0}\rangle$ obtained after convergence.

1.5 E.   Proof that the solution of the Ghost-GA equations is disentangled

from the auxiliary ghost subsystem

As explained in the main text, the Ghost-GA theory is formulated by embedding the Hamiltonian $\hat{H}$ within an extended Hilbert space, which includes also auxiliary (ghost) Fermionic degrees of freedom. Such a procedure is obviously licit, as the inequality underlying the variational principle, i.e.,

[TABLE]

where $E_{0}$ is the ground state energy of $\hat{H}$, remains valid even if $|\Psi\rangle$ is a generic state of the extended Hilbert space.

Since $\hat{H}$ lives only within the physical subsystem, we already know a-priori that the exact ground state (which realizes the above variational minimum) has to be disentangled from the auxiliary degrees of freedom. Here we demonstrate that, indeed, this condition is satisfied exactly by our converged Ghost-GA solution $|\Psi_{G}\rangle=\hat{\mathcal{P}}_{G}|\Psi_{0}\rangle$.

Let us assume that among the $\tilde{\nu}$ Fermionic degrees of freedom in $\hat{H}$, see Eq. (1), only $\nu$ are physical, while the other $\tilde{\nu}-\nu$ degrees of freedom are ghost modes. In other words, $\hat{H}$ is constructed utilizing only the physical ladder operators $\{c^{\dagger}_{R\alpha},c^{\phantom{\dagger}}_{R\alpha}\,|\,\alpha\in\{1,...,\nu\}\}$ (which implies, in particular, that $[\epsilon_{k}]_{\alpha\beta}=0\,\forall\,k,\forall\,\alpha,\beta>\nu$).

It is convenient to represent the Gutzwiller local maps $\hat{\mathcal{P}}_{R}$ as follows:

[TABLE]

where

[TABLE]

As we are going to show, since $\hat{H}$ depends only on the physical degrees of freedom, i.e., it is constructed utilizing only the physical operators $\{c^{\dagger}_{R\alpha},c^{\phantom{\dagger}}_{R\alpha}\,|\,\alpha\in\{1,...,\nu\}\}$, the coefficients defining the Gutzwiller projector are of the form:

[TABLE]

Consequently, the converged Ghost-GA solution $|\Psi_{G}\rangle=\hat{\mathcal{P}}_{G}|\Psi_{0}\rangle$ can be represented as

[TABLE]

where $|\Psi^{phys}_{G}\rangle$, which resides entirely within the physical subsystem, and $|\Psi^{ghost}_{G}\rangle$, which resides entirely within the ghost subsystem, are disentangled.

Let us now demonstrate Eq. (29). Since we assumed that $[\epsilon_{k}]_{\alpha\beta}=0\,\forall\,k,\forall\,\alpha,\beta>\nu$, from Eq. (20) it follows that

[TABLE]

Consequently, the ground state $|\Phi\rangle$ of $\hat{H}_{\text{emb}}$, see Eq. (11), is such that the subsystem generated by the ghost degrees of freedom $\{c^{\dagger}_{R\alpha},c^{\phantom{\dagger}}_{R\alpha}\,|\,\alpha>\nu\}$ is disentangled from the rest of the embedding system. Thus, from Eq. (12) it follows that the matrix $\tilde{\phi}$ can be represented as

[TABLE]

The proof of Eq. (29) (and, in turn, of Eq. (30)) follows immediately from Eqs. (32) and (13).

In summary, we have shown that the Ghost-GA solution, see Eq. (2), does not have any spurious entanglement with the auxiliary ghost degrees of freedom (as expected).

1.6 F.   Exploiting Eq. (30) numerically

Since we know a-priori that Eq. (30) must be satisfied by the converged result, it is computationally convenient to impose from the onset this condition onto the Ghost-GA variational parameters.

• •

The first simplification due to Eq. (30) arises from the fact that from Eq. (31) implies that

[TABLE]

Consequently, from Eq. (23) it follows that

[TABLE]

This equation enables us to reduce substantially the computational complexity of the root problem [Eq. (25)], as the dimension of the space of matrices satisfying Eq. (34) has only dimension $\nu\tilde{\nu}$ instead of $\tilde{\nu}^{2}$, see Eq. (18).

• •

The second simplification arises from the fact that, because of Eq. (31), the embedding Hamiltonian [Eq. (11)] is such that its subsystem generated by $\{\hat{c}^{\dagger}_{R\alpha},\hat{c}^{\phantom{\dagger}}_{R\alpha}\,|\,\alpha\in\{1,...,\nu\}\}\cup\{\hat{f}^{\dagger}_{Ra},\hat{f}^{\phantom{\dagger}}_{Ra}\,|\,a\in\{1,...,\tilde{\nu}\}\}$ is disentangled from the ghost degrees of freedom $\{\hat{c}^{\dagger}_{R\alpha},\hat{c}^{\phantom{\dagger}}_{R\alpha}\,|\,\alpha>\nu\}$. Consequently, it can be solved independently.

These observations reduce exponentially the scaling of the computational complexity of the problem as a function of the number of ghost orbitals used in the calculation.

1.7 G.   Similarities with theory of Matrix Product States (MPS) and

Projected Entangled Pair States (PEPS)

In the previous subsections we have shown that the Ghost-GA wavefunction $|\Psi_{G}\rangle=\hat{\mathcal{P}}_{G}|\Psi_{0}\rangle$ is disentangled from the auxiliary subsystem. However, as we pointed out in the main text, $|\Psi_{0}\rangle$ resides within the entire Hilbert space — including the ghost degrees of freedom. In other words, the Gutzwiller operator $\hat{\mathcal{P}}_{G}=\prod_{R}\hat{\mathcal{P}}_{R}$ maps $|\Psi_{0}\rangle$ — which is entangled with the auxiliary subsystem — into the a disentangled state represented as in Eq. (30). It is interesting to note that, in this respect, the Ghost-GA variational framework is very similar to the theory of Matrix Product States (MPS) and Projected Entangled Pair States (PEPS) [Verstraete et al., 2004; Schollwock, 2011].

As we are going to show, the main differences between the Ghost-GA theory and MPS/PEPS are the following:

• •

Within the Ghost-GA theory, the local maps (i.e., the Gutzwiller projectors $\hat{\mathcal{P}}_{R}$) act on a generic ”virtual” Slater determinant $|\Psi_{0}\rangle$ belonging to an extended Hilbert space, which is determined variationally.

• •

Instead, MPS/PEPS are variational techniques formulated by applying local maps to fixed ”pair states”, which are also constructed within an extended Hilbert space [Verstraete et al., 2004].

In this respect, the variational ansatz represented in Eq. (2) constitutes an extension of the MPS/PEPS variational space. Of course, as we pointed out above and in the main text, in our work we have also assumed the Gutzwiller Approximation [Gutzwiller, 1965] (which is exact only in the limit of infinite dimensions) and the Gutziller constraints, see Eq. (8) — which are approximations not employed in MPS/PEPS theory. Because of this reason, while our approach constitutes a variationally improved scheme with respect to the ordinary GA, it is not ”numerically exact”.

In order to illustrate more clearly the connection between the Ghost-GA theory and MPS/PEPS, let us consider a generic 1-dimensional Fermionic system belonging to a Hilbert space represented as $\mathcal{H}:=\bigotimes_{R=1}^{\Omega}\mathcal{H}_{R}$, where each local subsystem $\mathcal{H}_{R}$ is generated by a set of Fermionic operators $\{f^{\dagger}_{R\alpha}\}$, where $a\in\{1,...,{\tilde{\nu}}\}$. In this subsection we are going to assume also that $\tilde{\nu}$ is even. Following, e.g., Ref. [Verstraete et al., 2004], let us consider the following state:

[TABLE]

which is a tensor product of maximally entangled ”bonds” represented as

[TABLE]

where $f^{\dagger}_{R\,\Omega+1}=f^{\dagger}_{R1}$. Note that the state $|\Psi_{0}\rangle$ defined in Eq. (35) is a Slater determinant.

Let us now assume that only $\nu<\frac{\tilde{\nu}}{2}$ of the degrees of freedom are physical, i.e., that we are interested in the ground state of a Hamiltonian $\hat{H}$ constructed only with the operators $\{c^{\dagger}_{R\alpha},c^{\phantom{\dagger}}_{R\alpha}\,|\,\alpha\in\{1,...,\nu\}\}$, which are connected to the operators $\{f^{\dagger}_{Ra},f^{\phantom{\dagger}}_{Ra}\,|\,a\in\{1,...,\tilde{\nu}\}\}$ by unitary transformations

[TABLE]

As in Sec. I E, in order to construct a variational approximation of the ground state of $\hat{H}$, we apply on $|\Psi_{0}\rangle$ a Gutzwiller operator $\hat{\mathcal{P}}_{G}=\prod_{R}\hat{\mathcal{P}}_{R}$, see Fig. 1, where $\hat{\mathcal{P}}_{R}$ is the most general map acting over all of the local degrees of freedom labeled by $R$ represented as follows:

[TABLE]

where

[TABLE]

and the coefficients defining the Gutzwiller projector are of the form:

[TABLE]

Following Ref. [Verstraete et al., 2004], it can be readily verified by inspection that the so obtained state $|\Psi_{G}\rangle=\hat{\mathcal{P}}_{G}|\Psi_{0}\rangle$ can be represented as follows:

[TABLE]

where $|\Psi^{phys}_{G}\rangle$ is a MPS with bond dimension $2^{\tilde{\nu}}$, which resides entirely within the physical subsystem.

In particular, the observation above clarifies that, within the Ghost-GA variational ansatz, increasing the number of ghost orbitals amounts essentially to increase the so-called ”bond dimension” [Verstraete et al., 2004], which is the parameter controlling the accuracy of the variational ansatz also in MPS/PEPS. However, we remark that within the Ghost-GA theory the state $|\Psi_{0}\rangle$ is variationally determined, i.e., it is not restricted to the form [Eq. (36)]. In future works, it would be interesting to apply the ghost-orbital Hilbert space extension without employing the Gutzwiller constraints and the Gutzwiller approximations, e.g., utilizing the Variational Monte Carlo method [Ceperley et al., 1977].

2 II.   Ghost-GA excitations from variational parameters

2.1 A.   Gutzwiller excitations and ARPES spectra in the standard multi-orbital GA

Let us assume that the Gutzwiller energy minimum is realized by a solution of Eqs. (19)-(24) identified by the variational parameters $|\Psi_{0}\rangle,E,|\Phi\rangle,E^{c},\mathcal{R},\lambda,\mathcal{D},\lambda^{c},\Delta$. Thus, $|\Psi_{0}\rangle$ is the ground state of $\hat{H}_{\text{qp}}$ corresponding to the parameters $\mathcal{R},\lambda$, see Eq. (10). Such a solution corresponds to a Gutzwiller wavefunction which can be represented as

[TABLE]

see Eq. (2). It is important to observe that, in the thermodynamical limit $\mathcal{N}\rightarrow\infty$, Eqs. (19)-(24) are satisfied also by $\psi^{\dagger}_{kn}|\Psi_{0}\rangle,E+{\epsilon}^{*}_{kn},|\Phi\rangle,E^{c},\mathcal{R},\lambda,\mathcal{D},\lambda^{c},\Delta$, where $\psi^{\dagger}_{kn}$ is any eigenoperator of $\hat{H}_{\text{qp}}[\mathcal{R},\lambda]$, see Eq. (10), and ${\epsilon}^{*}_{kn}$ is the corresponding eigenvalue. Each one of the corresponding Gutzwiller states can be represented as

[TABLE]

Similarly, Eqs. (19)-(24) are satisfied also by $\psi^{\phantom{\dagger}}_{kn}|\Psi_{0}\rangle,E-{\epsilon}^{*}_{kn},|\Phi\rangle,E^{c},\mathcal{R},\lambda,\mathcal{D},\lambda^{c},\Delta$, which correspond to Gutzwiller states represented as

[TABLE]

Since the states [Eqs. (43) and (44)] correspond to saddle points of the Gutzwiller energy functional, they can be interpreted as variational approximations to many-body excitations of the system (with a number of particles differing by 1 with respect to the ground state). It is important to note that each one of these states is connected to a single-particle excitation of $\hat{H}_{\text{qp}}$ through the action of the many-body operator $\hat{\mathcal{P}}_{G}$. Thus, these states represent complex collective excitations of the system, which are commonly called ”Gutzwiller-Landau quasiparticles” [Bünemann et al., 2003].

As pointed out in previous publications, see, e.g., Ref. [Lanatà, 2009] and the references therein, this information allows us to evaluate the Gutzwiller ARPES spectra, which is defined as follows:

[TABLE]

The main idea in order to evaluate approximately Eq. (45) consists in inserting projectors over the Gutzwiller-Landau quasiparticle states $\sum_{n}|\Psi^{p}_{Gkn}\rangle\langle\Psi^{p}_{Gkn}|$ and $\sum_{n}|\Psi^{h}_{Gkn}\rangle\langle\Psi^{h}_{Gkn}|$ (see Eqs. (43) and (44)) in Eq. (45) and use the following identities:

[TABLE]

where $f^{\phantom{\dagger}}_{ka}$ are the operators appearing in Eq. (10). From the identities above one can readily obtain [Lanatà, 2009]:

[TABLE]

Note that, in general, [Eqs (44) and (44)] are not a complete basis. In fact, we have that, within the approximations employed in Eq. (48),

[TABLE]

Within the standard GA theory, the interpretation of $Z$ as the matrix of quasiparticle weights is motivated by the fact that the Green’s function corresponding to Eq. (48) is

[TABLE]

Thus, if $\mathcal{R}$ is invertible, it can be straightforwardly verified that:

[TABLE]

It is very important to note that Eq. (51), which provides a linear expression for the self-energy as a function of the frequency $\omega$, is valid only if $\mathcal{R}$ is invertible. However, we already know that this condition is not satisfied within the Ghost-GA theory, see Eq. (34). This important point will be further stressed in the next subsection.

2.2 B.   Gutzwiller excitations and ARPES spectra in the Ghost-GA theory

As explained in the main text and in the previous section of the present supplemental material, the Ghost-GA theory consists in applying the ordinary multi-orbital Gutzwiller theory to the Hamiltonian ”embedded” within the extended Hilbert space obtained by introducing auxiliary ghost orbitals. The only difference with respect to the ordinary theory is that, as shown in the main text, utilizing the enlarged Gutzwiller variational space leads to a better approximation to the ground state $|\Psi_{G}\rangle$ and to the excited states $|\Psi^{p}_{Gkn}\rangle$ and $|\Psi^{h}_{Gkn}\rangle$, see Eqs. (43) and (44), which are all saddle points of the Ghost-GA energy functional.

Let us consider again the Hamiltonian $\hat{H}$, see Eq. (1), embedded within an enlarged Hilbert space. In other words, we assume that $\hat{H}$ is constructed utilizing only the physical ladder operators $\{c^{\dagger}_{R\alpha},c^{\phantom{\dagger}}_{R\alpha}\,|\,\alpha\in\{1,...,\nu\}\}$ (which implies, in particular, that $[\epsilon_{k}]_{\alpha\beta}=0\,\forall\,k,\forall\,\alpha,\beta>\nu$).

The Ghost-GA spectral function of the system is given by the components of the spectral function

[TABLE]

such that $\alpha,\beta\in\{1,...,\nu\}$. The same steps utilized in the previous section lead to the following approximation to the physical Green’s function:

[TABLE]

Note also that, as expected, from Eqs. (34) and (53) it follows that

[TABLE]

i.e., only the physical components of the Green’s function are non-zero.

The Ghost-GA physical self-energy is given, by definition, by the Dyson equation, i.e., it is the function $\Sigma(\omega)$ satisfying the following equation:

[TABLE]

where $G_{\alpha\beta}(\omega)$ is given by Eq. (53). It is important to note that, as pointed out above and in the main text, $\mathcal{R}$ is not invertible, see Eq. (34). Consequently, the Ghost-GA self-energy $\Sigma(\omega)$ defined above is not necessarily a linear function of the frequency $\omega$.

Remarkably, Eq. (55) provides us with an analytical expression for the Ghost-GA spectral function $A(k,\omega)$ in terms of the bare physical dispersion $\epsilon_{k}$ and a few renormalization parameters, i.e., the matrices $\tilde{\mathcal{R}}$ and $\tilde{\lambda}$. Of course, the entries of $\tilde{\mathcal{R}}$ and $\tilde{\lambda}$ depend on the specific system considered, as they have to be calculated numerically from Eqs (19)-(24). In particular, Eqs. 8-13 of the main text have been obtained from Eq. (55).

Note that, within the context of the Ghost-GA theory, the excited states $|\Psi^{p}_{Gkn}\rangle$ and $|\Psi^{h}_{Gkn}\rangle$ are constructed by applying the Gutzwiller projector to states belonging to the extended Hilbert space. In this respect, these mathematical form of these states present suggestive formal analogies with the ”hidden Fermi liquid” excitations previously introduced by P. W. Anderson [Anderson, 2008; Jain and Anderson, 2009].

We point out also that the construction above is not specific to the single band Hubbard Hamiltonian considered in this work, as it can be straightforwardly applied to generic Hubbard Hamiltonian — with an arbitrary number of physical orbitals.

3 III.   Additional Remarks about the benchmark calculations of the single band Hubbard model

3.1 A.   Behavior of the Ghost-GA self-energy

As discussed in the main text, the Ghost-GA self-energy of the single band Hubbard Hamiltonian can be represented as follows in the metallic phase:

[TABLE]

while it can be represented as follows in the Mott phase:

[TABLE]

The values of the parameters $z$ and $h$ evolve as illustrated in Fig. 2 of the main text. The behavior of the Ghost-GA self-energy is shown in Fig. 2 in comparison with DMFT.

As discussed in the main text, some of the main features of the self energy are captured by Eqs. (56) and (57). The differences with respect to the DMFT solution arise mostly from the fact that within the Ghost-GA theory, by construction, the self-energy can develop only poles, see Eqs. (56) and (57), but can not capture branch-cut singularities on the real axis. In other words, the so-called ”scattering rate” is not captured by our approximation (as it is not captured by the ordinary GA).

It is interesting to observe that in the limit of $U\rightarrow 0$ we have $z\rightarrow 1$ and $h\rightarrow 0$. Thus, from Eq. (56) we deduce that $\lim_{U\rightarrow 0}\Sigma(\omega)\rightarrow 0\;\forall\,\omega$, i.e., the uncorrelated limit is described exactly by our theory. Note that this result is to be expected as, by construction, the Ghost-GA variational space includes all of the Slater determinants.

3.2 B.   Ghost-GA solution utilizing 3 ghost orbitals

In this work, the accuracy of the Ghost-GA theory has been explicitly demonstrated (in all regimes of interaction strength) on the single-band Hubbard Model for a semi-circular density of states (DOS) — which corresponds, e.g., to the Bethe lattice in the limit of infinite coordination number, where DMFT is exact. Specifically, the calculations discussed in the main text were performed using 2 ghost orbitals in the metallic phase and 1 ghost orbital in the Mott phase, as this was the minimal ghost extension enabling us to capture the main features of the ARPES spectra and of the $T=0$ phase diagram of the system. In this section we are going to discuss how the Ghost-GA result depends on the number of ghost orbitals used in the calculations.

The Ghost-GA physical Green’s function of the single-band Hubbard model

[TABLE]

see the third member of Eq. (53), is analytical by construction everywhere over the upper-half plane, and has poles over the eigenvalues $\tilde{\epsilon}^{*}_{kn}$ of:

[TABLE]

However, the physical spectral weights of these poles is generally not $1$ because of the presence of the renormalization coefficients $\tilde{\mathcal{R}}$ in Eq. (58).

As mentioned in the main text, if 2 ghost orbitals are used to study the Mott phase of the Hubbard model, all of the results, including both the ground state properties and the ARPES spectra, remain essentially unchanged. In particular, even though $\hat{H}_{\text{qp}}$ has 3 bands (by construction), the converged values of the renormalization parameters are such that one of these bands is flat and acquires spectral weight [math] (i.e., it does not contribute to the ARPES spectra), while the other 2 bands do not display any visible difference with respect to those computed using only 1 ghost orbital. Here we show that, as expected, the same mechanism occurs when 3 ghost orbitals are used to describe the system, both in the metallic and in the insulating phases of the single-band Hubbard model with semicircular DOS at half-filling.

In Figs. 4 and 4 the Ghost-GA energy-resolved spectral function $A(\epsilon,\omega)=-\frac{1}{\pi}\text{Im}\,{G(\epsilon,\omega)}$, see Eq. (58) is analyzed for 2 values of $U$ in comparison with the Ghost-GA quasiparticle spectral function \mathrm{Tr}\big{[}\delta(\omega-[\tilde{\mathcal{R}}\,\tilde{\epsilon}\,\tilde{\mathcal{R}}^{\dagger}+\tilde{\lambda}])\big{]}=\sum_{n}\delta(\omega-\tilde{\epsilon}^{*}_{n}), see Eq. (59). The Ghost-GA results obtained in the main text — utilizing 2 ghost orbitals in the metallic phase and 1 ghost orbital in the Mott phase — are compared with the results obtained using 3 ghost orbitals for both phases.

This analysis shows that, in all cases analyzed so far, increasing the number of orbitals (ghost plus physical) beyond the minimal number necessary in order to match the number of bands in the system does not affect the physical spectral function $\mathcal{A}(\epsilon,\omega)$. In particular, in the results shown in Figs. 4 and 4 we observe that the physical spectral function obtained with 3 ghost orbitals displays 1 flat band with zero spectral weight in the metallic phase, while it displays 2 degenerate flat bands with 0 spectral weight in the Mott phase.

3.3 C.   Benchmark calculations of the single-band Hubbard model away from half-filling

In order to study the Hubbard Hamiltonian away from half filling, here we work in the grand-canonical ensemble, i.e., we study the following Hamiltonian:

[TABLE]

Note that Eq. (60) is represented in such a way that the particle-hole symmetric case (half-filling) corresponds to $\mu=0$.

In Fig. 5 is shown the evolution as a function of the chemical potential $\mu$ of the total energy, the local double occupancy, the QP weight $z$ and the local occupancy, for 2 values of the Hubbard interaction strength $U$. The Ghost-GA results are shown in comparison with the ordinary GA theory and with DMFT.

As in the case of half filling considered above and in the main text, the agreement between Ghost-GA and DMFT is quantitatively remarkable.

4 IV.   Ghost-GA benchmark calculations of the Single-orbital Anderson Impurity Model

In order to provide further evidence of the quality the Ghost-GA theory and, in turn, to further validate the results provided in the main text, in this section we provide also benchmark calculations of the single-band Anderson Impurity Model (AIM).

[TABLE]

where $\sigma$ is the spin label. In particular, we are going to consider the case of a flat density of states, whose hybridization function [Hewson, 1997] can be represented as follows:

[TABLE]

where $\theta$ is the Heaviside step function.

In Fig. 6 is shown the evolution as a function of the Hubbard interaction strength $U$ of the impurity double occupancy and the QP weight $z$. The Ghost-GA results obtained by adding 2 ghost orbitals on the impurity site are shown in comparison with the ordinary GA theory and with NRG.

Consistently with the benchmark calculations of the Hubbard model considered above and in the main text, here we observe that the accuracy of the Ghost-GA solution is quantitatively remarkable also for the AIM. In particular, we note that the Ghost-GA results are substantially more accurate with respect to the ordinary GA — as expected, since the Ghost-GA theory enlarges the ordinary Gutzwiller variational space.

It is interesting to analyze more in detail the behavior of $z$ (which is commensurate with the Kondo temperature [Hewson, 1997]), in the strong-coupling regime $U\gg\Gamma$. It is well known that the ordinary Gutzwiller theory provides, even without ghost orbitals, an exponentially decaying Kondo scale in the large-$U$ limit: $z\sim e^{-\frac{\pi}{16}\frac{U}{\Gamma}}$. The only well-known discrepancy of this result with respect to the exact solution of the AIM is that the Gutzwiller prefactor at the exponent is $\pi/16$, while the exact universal numerical value is $\pi/8$, see, e.g., Eqs. 101-102 of Ref. [Lanatà, 2010] and the references therein. A numerical analysis of the results displayed in Fig. 6 shows that the Ghost-GA theory obtained by adding 2 ghost orbitals provides, in particular, a more accurate estimation of the universal prefactor in front of $U/\Gamma$ at the exponent of the expression for $T_{K}$, which behaves as $z\sim e^{-\frac{\pi}{13.7}\frac{U}{\Gamma}}$ in the large-$U$ regime.

A more detailed study of the AIM will be provided in future publications.

5 V.   Ghost Rotationally Invariant Slave Boson (Ghost-RISB) theory

As mentioned in the conclusions of the main text, the Ghost-GA theory can be equivalently reformulated as the mean field approximation of an alternative formulation of the Rotationally Invariant Slave Boson (RISB) theory, which we call Ghost-RISB theory. For completeness, this construction is briefly described in this section.

Let us consider a generic Hubbard Hamiltonian represented as in Eq. (1). In Ref. [Lanatà et al., 2017] it was demonstrated that Eq. (1) can be equivalently reformulated as a gauge theory in an auxiliary Hilbert space (i.e., the RISB theory), and that this exact reformulation of the many-body problem reduces to the ordinary GA at the mean-field level.

Here we observe that the mapping mentioned above could be equivalently applied to the Hubbard Hamiltonian expressed within a Ghost-GA extended Hilbert space. Of course, this construction would result in an alternative exact reformulation of the many body problem, which we call Ghost-RISB theory. However, by construction, the Ghost-RISB theory reduces to the Ghost-GA theory at the mean field level — which is substantially more accurate with respect to the ordinary GA. Because of this reason, it would be very interesting to explore the possibility of taking into account fluctuations beyond mean-field within the framework of the Ghost-RISB, as it has been done in previous works within the ordinary SB theory [Dao and Frésard, 2017].