Photoemission Spectra from the Extended Koopman’s Theorem, Revisited

1 Introduction

The Extended Koopman’s Theorem (EKT) (Morrell et al., 1975; Smith and Day, 1975) has been derived in quantum chemistry and used within various frameworks, from functional theories based on reduced quantities, such as reduced-density matrix functional theory (Gilbert, 1975) (e.g., Pernal and Cioslowski, 2005; Leiva and Piris, 2005; Piris et al., 2012; Piris et al., 2013) and many-body perturbation theory based on Green’s functions (Hedin, 1965) (e.g., Dahlen and van Leeuwen, 2005; Stan et al., 2006; Stan et al., 2009), to wavefunction-based methods (e.g., Cioslowski et al., 1997; Kent et al., 1998; Bozkaya, 2013; Zheng, 2016; Bozkaya and Ünal, 2018; Pavlyukh, 2019; Lee et al., 2021). The EKT allows one to calculate energies corresponding to charged excitations. Although it can be formulated both for ionization potentials (IPs) and electron affinities (EAs), it has been widely used only for the former, whereas for the latter applications have been limited to the calculation of the lowest EA as the first IP of the (N+1)-electron system (in case of finite systems), with N the number of electrons in the reference system. There exist hence many benchmarks for the IPs. So far, the method has been mainly used for finite systems. The EKT is known to be in principle (i.e., using exact ingredients, namely the one- and two body density matrices, as we shall see) exact for the first ionization potential (Katriel and Davidson, 1980; Sundholm and Olsen, 1993). In the solid state instead there are only a few applications which only focus on the band structure. It would be desirable to have also the spectral function, i.e., the spectrum of electron addition and removal energies weighted by the Dyson amplitudes, which measure the overlap between the eigenstates of the (N + 1) − electron ((N − 1) − electron) system and the ground state of the N-electron system where an electron has been added (removed). The spectral function is related to photoemission spectroscopy, which gives precious information about the electronic structure and excitations in a system, and, moreover, allows one to study metal-insulator transitions, of paramount importance in condensed matter.

A simple way to calculate the spectral function A(ω) is through the imaginary part of the one-body Green’s function G(ω), as $A(\omega )=\frac{1}{\pi }\text{sgn}(\mu -\omega )\mathfrak{I}G(\omega )$, where μ is the chemical potential. The one-body Green’s function is the fundamental quantity of many-body perturbation theory; it can be obtained from the Dyson equation G = G₀ + G₀ΣG, in which G₀ and G are the noninteracting and interacting Green’s functions, respectively, and Σ the so-called self-energy, which contains all the many-body effects of the system. This latter quantity needs to be approximated in practical calculations. Commonly used approximations, such as the well-known GW approximation (Hedin, 1965), cannot capture the Mott physics (Romaniello et al., 2009; Romaniello et al., 2012; Di Sabatino et al., 2015; Di Sabatino et al., 2016; Di Sabatino et al., 2021). Therefore much effort is devoted to develop better approximations to Σ (Springer et al., 1998; Zhukov et al., 2004; Shishkin et al., 2007; Kuneš et al., 2007; Guzzo et al., 2011; Romaniello et al., 2012; Lischner et al., 2013; Stefanucci et al., 2014) or to develop novel ways to determine G (Lani et al., 2012; Berger et al., 2014). In this spirit in these last years we have developed the many-body effective energy theory (MEET) (Di Sabatino et al., 2016), in which the spectral function is expressed in terms of density matrices, or, alternatively, in terms of moments of G, as reported in Ref. (Di Sabatino et al., 2019). This has allowed us to describe the band gap in several paramagnetic transition-metal oxides (Di Sabatino et al., 2016; Di Sabatino et al., 2019; Di Sabatino et al., 2021), such as NiO, which are considered strongly correlated materials and which are described as metals by static mean-field theories, such as DFT, and by GW. This is an important result. However the band gap is hugely overestimated by the MEET within the current low-order approximation in terms of the (approximate) one- and two-body density matrices. Improvements are needed, either by going to higher-order density matrices, which, however, is not guaranteed to converge, or by introducing some sort of screening in the equations. Recently we have obtained promising results for the description of the insulator-to-metal transition of PM FeO under pressure by combining the MEET and the local-density approximation (LDA) (Di Sabatino et al., 2021), and we are currently working on introducing electron-hole screening in the MEET equations. However there is another path which we can explore, and this comes from the relation between the MEET and the EKT. As we will show in the following, within a given basis, the removal and addition energies obtained within the MEET at the lowest-order approximation are equal to the EKT removal and addition energies within the diagonal approximation. In this work we discuss this link and its impact on both theories.

The paper is organized as follows. In Theory we give the basic equations of the EKT and the MEET and we make the link between them. The Hubbard dimer and bulk silicon are used to illustrate the difference between the EKT and the MEET (removal and addition) energies in Mutual Insights and Illustration. In Conclusions and Perspectives we draw our conclusions and perspectives.

2 Theory

In this section we briefly review the MEET and EKT methods, and we make the link between the two. We will consider an N-electron system governed by the following Hamiltonian in second quantization

$\widehat{H}=\sum _{ij}{h}_{ij}{\widehat{a}}_i^{†}{\widehat{a}}_j+\frac{1}{2}\sum _{ijkl}{V}_{ijkl}{\widehat{a}}_i^{†}{\widehat{a}}_j^{†}{\widehat{a}}_l{\widehat{a}}_k,$

where $\widehat{a}$ and ${\widehat{a}}^{†}$ are the annihilation and creation operator, respectively, $h_{ij}=\int d\mathbf{x}{\varphi }_i^{*}(\mathbf{x})h(\mathbf{r}){\varphi }_j(\mathbf{x})$ are the matrix elements of the one-particle noninteracting Hamiltonian h(r) = − ∇²/2 + v_ext(r), with v_ext an external potential, and $V_{ijkl}=\int d\mathbf{x}d{\mathbf{x}}^{\prime }{\varphi }_i^{*}(\mathbf{x}){\varphi }_j^{*}({\mathbf{x}}^{\prime })v_c(\mathbf{r},{\mathbf{r}}^{\prime }){\varphi }_k(\mathbf{x}){\varphi }_l({\mathbf{x}}^{\prime })$ are the matrix elements of the Coulomb interaction v_c. Here x = (r, α) combines space and spin variables and i, j, … denote both space and spin labels (they will be made explicit only when necessary).

2.1 Key Equations of the MEET

Within the MEET the time-ordered 1-body Green’s function G(ω) at zero temperature is split into removal (R) and addition (A) parts as G(ω) = G^R(ω) + G^A(ω). In the following we concentrate on the diagonal elements of G, which are related to photoemission spectra. Within the MEET the diagonal matrix elements of G^R/A(ω) are written in terms of an effective energy ${\delta }_i^{R/A}(\omega )$ as (Di Sabatino et al., 2016):

$\begin{array}{ccc}\hfill G_{ii}^R(\omega )& =\hfill & \frac{{\gamma }_{ii}}{\omega -{\delta }_i^R(\omega )-\text{i}\eta },\hfill \end{array}(1)$

$\begin{array}{ccc}\hfill G_{ii}^A(\omega )& =\hfill & \frac{1-{\gamma }_{ii}}{\omega -{\delta }_i^A(\omega )+\text{i}\eta },\hfill \end{array}(2)$

with γ_ii the diagonal matrix element of the one-body density matrix in a given basis set. We note that a similar effective energy can be introduced also for the off-diagonal elements of G^R/A. The spectral function is hence expressed as

$A_{ii}(\omega )={\gamma }_{ii}\delta (\omega -{\delta }_i^R(\omega ))+(1-{\gamma }_{ii})\delta (\omega -{\delta }_i^A(\omega )),(3)$

where the symbol δ on the right-hand side indicates the Dirac delta function. In our previous works we have chosen the basis set of natural orbitals, i.e., the orbitals which diagonalize the one-body reduced density matrix. In this case γ_ii = n_i, i.e., the natural occupation numbers. This choice has been made based on our results on exactly solvable Hubbard clusters, where the MEET performs very well. (Di Sabatino et al., 2016) However this does not guarantee that it remains the best choice for more realistic systems. In fact this is not the case as we shall see.

The effective energy ${\delta }_i^{R/A}(\omega )$ can be written as an expansion in terms of reduced density matrices. The expression truncated at the level of the one- and two-body reduced density matrices (2-RDM) reads (in the basis of natural orbitals).

$\begin{array}{ccc}\hfill {\delta }_i^{R,(1)}& =\hfill & h_{ii}+\frac{1}{n_i}\sum _{klm}{V}_{imkl}{\mathrm{\Gamma }}_{klmi}^{(2)}\hfill \end{array}(4)$

$\begin{array}{ccc}\hfill {\delta }_i^{A,(1)}& =\hfill & h_{ii}+\frac{1}{(1-n_i)}\sum _k\left(V_{ikik}-V_{ikki}\right)n_k\hfill \\ \hfill & \hfill & -\frac{1}{(1-n_i)}\sum _{klm}{V}_{imkl}{\mathrm{\Gamma }}_{klmi}^{(2)},\hfill \end{array}(5)$

where ${\mathrm{\Gamma }}_{klmi}^{(2)}=⟨{\mathrm{\Psi }}_0^N|{\widehat{a}}_i^{†}{\widehat{a}}_m^{†}{\widehat{a}}_l{\widehat{a}}_k|{\mathrm{\Psi }}_0^N⟩$ are the matrix elements of the two-body reduced density matrix, with ${\mathrm{\Psi }}_0^N$ the ground-state wavefunction of the N-electron system. As discussed in Ref. (Di Sabatino et al., 2016) the various approximations ${\delta }_i^{R/A,(n)}(\omega )$ are related to the n-th moments

${\mu }_{n,i}^{R/A}=\frac{{\sum }_k{B}_{ii}^{k,R/A}{({ϵ}_k^{R/A})}^n}{{\sum }_k{B}_{ii}^{k,R/A}}$

of the $G_{ii}^{R/A}(\omega )$. Here ${ϵ}_k^R=(E_0^N-E_k^{N-1})$ and ${ϵ}_k^A=(E_k^{N+1}-E_0^N)$ are removal and addition energies, respectively, and

$\begin{array}{ll}\hfill B_{ii}^{k,R}& =\langle {\mathrm{\Psi }}_0^N|{\widehat{c}}_i^{†}|{\mathrm{\Psi }}_k^{N-1}\rangle \langle {\mathrm{\Psi }}_k^{N-1}|{\widehat{c}}_i|{\mathrm{\Psi }}_0^N\rangle \hfill \\ \hfill B_{ii}^{k,A}& =\langle {\mathrm{\Psi }}_0^N|{\widehat{c}}_i|{\mathrm{\Psi }}_k^{N+1}\rangle \langle {\mathrm{\Psi }}_k^{N+1}|{\widehat{c}}_i^{†}|{\mathrm{\Psi }}_0^N\rangle ,\hfill \end{array}$

with $E_0^N$ and ${\mathrm{\Psi }}_0^N$ the ground-state energy and wave function of the N-electron system and $E_k^{N\pm 1}$ and ${\mathrm{\Psi }}_k^{N\pm 1}$ the kth state energy and wave function of the (N ± 1)-electron system. This allows for a more compact expression of $G_{ii}^{R/A}(\omega )$ as a continued fraction of moments

$G_{ii}^R=\frac{n_i}{\omega -{\mu }_{1,i}^R\frac{\omega -{\mu }_{1,i}^R\dots }{\omega -\frac{{\mu }_{2,i}^R}{{\mu }_{1,i}^R}\dots }},(6)$

(and similarly for $G_{ii}^A$). More details on the continued fraction expression for G can be found in Refs (Di Sabatino et al., 2016; Di Sabatino et al., 2019). At the level of δ^R/A,(1), the Green’s function depends only on the first moment, while neglecting all the higher-order frequency-dependent corrections. As shown in Ref. (Di Sabatino et al., 2019) this means that each component $G_{ii}^{R/A}$ has only one pole which is a weighted average of all the poles of $G_{ii}^{R/A}$. If each component of G has a predominant quasiparticle peak, this is a good approximation, provided that the approximation to the first moment is accurate enough. At the level of δ^R/A,(2) the Green’s function depends on the first and second moments; since now the corrections are frequency-dependent more poles appear (namely, two removal and two addition poles for each component of G, which are visible if the corresponding weights are nonzero). This approximation tends to reproduce the two most dominant removal/addition peaks for each component of G. Higher-order moments will produce more poles; however, approximations become quickly uncontrolled (Di Sabatino, 2016), which can lead to unphysical results.

2.2 Key Equations of the EKT

Within the EKT one starts from the following approximation for the removal energy ${ϵ}_i^R$ (Kent et al., 1998)

${ϵ}_i^R=-\frac{\langle {\mathrm{\Psi }}_0^N|{\widehat{O}}_i^{†}[\widehat{H},{\widehat{O}}_i]|{\mathrm{\Psi }}_0^N\rangle }{\langle {\mathrm{\Psi }}_0^N|{\widehat{O}}_i^{†}{\widehat{O}}_i|{\mathrm{\Psi }}_0^N\rangle }(7)$

with ${\mathrm{\Psi }}_0^N$ the ground-state many-body wave function of the N-electron system, and ${\widehat{O}}_i={\sum }_k{C}_{ki}^R{\widehat{a}}_k$, ${\widehat{O}}_i^{†}={\sum }_k{C}_{ki}^{R*}{\widehat{a}}_k^{†}$, with $\{C_{ki}^R\}$ a set of coefficients to be determined. The stationary condition (with respect to the coefficients $C_{ki}^R$) for ${ϵ}_i^R$ leads to the secular equation

$({\mathbf{V}}^R-{ϵ}_i^R{\mathbf{S}}^R){\mathbf{C}}_i^R=0,(8)$

with $V_{ij}^R=-⟨{\mathrm{\Psi }}_0^N|{\widehat{a}}_j^{†}[\widehat{H},{\widehat{a}}_i]|{\mathrm{\Psi }}_0^N⟩$ and S^R the one-body density matrix $S_{ij}^R={\gamma }_{ij}=⟨{\mathrm{\Psi }}_0^N|{\widehat{a}}_j^{†}{\widehat{a}}_i|{\mathrm{\Psi }}_0^N⟩$. If one defines the matrix ${\mathbf{\Lambda }}^R={[{\mathbf{S}}^R]}^{-1}{\mathbf{V}}^R$ in the basis of natural orbitals, with $S_{ij}^R=n_i{\delta }_{ij}$ and works out the commutator in $V_{ij}^R$, one arrives at

${\mathrm{\Lambda }}_{ij}^R=\frac{1}{n_i}\left[n_i{h}_{ji}+\sum _{klm}{V}_{jmkl}{\mathrm{\Gamma }}_{klmi}^{(2)}\right].(9)$

The eigenvalues of Λ^R are the removal energies. (Morrell et al., 1975; Pernal and Cioslowski, 2005) By comparing to Eq. 4 it becomes clear that the diagonal element of Λ^R are the removal energy of the MEET within the low-order approximation. The diagonal element of Λ^R are also referred in literature as the energies of the EKT within the diagonal approximation (DEKT).

Similar equations hold for the addition energies. One can indeed define the addition energy ${ϵ}_i^A$ as

${ϵ}_i^A=\frac{\langle {\mathrm{\Psi }}_0^N|[\widehat{H},{\widehat{O}}_i]{\widehat{O}}_i^{†}|{\mathrm{\Psi }}_0^N\rangle }{\langle {\mathrm{\Psi }}_0^N|{\widehat{O}}_i{\widehat{O}}_i^{†}|{\mathrm{\Psi }}_0^N\rangle }(10)$

and in a similar way as for ${ϵ}_i^R$ we arrive at the eigenvalue equation

$({\mathbf{V}}^A-{ϵ}_i^A{\mathbf{S}}^A){\mathbf{C}}_i^A=0,(11)$

with $V_{ij}^A=⟨{\mathrm{\Psi }}_0^N|{\widehat{a}}_i[\widehat{H},{\widehat{a}}_j^{†}]|{\mathrm{\Psi }}_0^N⟩$ and S^A related to the one-body density matrix as $S_{ij}^A=1-{\gamma }_{ij}$. Similarly to the removal energy problem, using the basis of natural orbitals, one can work out the commutator in $V_{ij}^A$ and reformulate the problem in terms of the matrix ${\mathbf{\Lambda }}^A={[{\mathbf{S}}^A]}^{-1}{\mathbf{V}}^A$¹, which reads

$\begin{array}{ccc}\hfill {\mathrm{\Lambda }}_{ij}^A& =\hfill & \frac{1}{(1-n_i)}\times \hfill \\ \hfill & \hfill & \left[(1-n_i)h_{ji}+\sum _k\left(V_{jkik}-V_{jkki}\right)n_k-\sum _{klm}{V}_{jmkl}{\mathrm{\Gamma }}_{klmi}^{(2)}\right].\hfill \end{array}(12)$

Again, the diagonal elements of Λ^A are the MEET addition energies within the approximation given in Eq. 5.²

3 Mutual Insights and Illustration

Now that we have established the link between the EKT and the MEET we will study how these theories can benefit from mutual insight.

3.1 Hubbard Dimer

We use a modified version of the Hubbard dimer in which the on-site Coulomb interaction is different for the two sites. Its Hamiltonian is given by

$H=-t\sum _{\begin{array}{c}i,j=\mathrm{1,2}\\ i\ne j\end{array}}\sum _{\sigma }{\widehat{a}}_{i\sigma }^{†}{\widehat{a}}_{j\sigma }+U_1{\widehat{n}}_{1\uparrow }{\widehat{n}}_{1\downarrow }+U_2{\widehat{n}}_{2\uparrow }{\widehat{n}}_{2\downarrow },(13)$

where i, j run over the sites, ${\widehat{n}}_{i\sigma }={\widehat{a}}_{i\sigma }^{†}{\widehat{a}}_{i\sigma }$, U_i is the on-site interaction at site i, − t is the hopping kinetic energy (the site energy ϵ₀ has been set to zero). Contrary to the standard dimer with a unique on-site interaction, in the case of two different on-site interactions the Λ^R/A are not diagonal in the basis of natural orbitals. Therefore, this model allows us to study the effect of the diagonalization on the removal/addition energies in the diagonal approximation. The model can represent the case of a heteronuclear diatomic molecule in a minimal basis set in which the valence orbitals of the two atoms are of different nature, such as HCl or NiO, for example. We note that also using the asymmetric Hubbard dimer with two different site energies the EKT equations are not diagonal in the basis of natural orbitals, however the difference between EKT and DEKT energies is not significant.

3.1.1 Insights Into the EKT

Making the parallel with the MEET, one can readily define the EKT spectral function as.

$\begin{array}{ll}\hfill A_{ii}^R(\omega )& ={\gamma }_{ii}\delta (\omega -{ϵ}_i^{EKT,R}),\hfill \end{array}(14)$

$\begin{array}{ll}\hfill A_{ii}^A(\omega )& =(1-{\gamma }_{ii})\delta (\omega -{ϵ}_i^{EKT,A}),\hfill \end{array}(15)$

with γ_ii and 1 − γ_ii the diagonal matrix element of the one-body density matrix in the basis which diagonalizes ${\mathrm{\Lambda }}_{ij}^R$ and ${\mathrm{\Lambda }}_{ij}^A$, respectively (not necessarily the same for ${\mathrm{\Lambda }}_{ij}^R$ and ${\mathrm{\Lambda }}_{ij}^A$). We notice that the factor γ_ii should refer to a proper one-body density matrix, i.e., a one-body density matrix which fulfils the ensemble N-representability constraints. In our case this is guaranteed by the total energy minimization (which includes the constraint 0 ≤ n_i ≤ 1) in RDMFT. Moreover, as for the MEET removal (addition) energies (in its lowest-order approximation), the removal (addition) EKT energies can be interpreted in terms of the first moment of $G_{ii}^R$ ($G_{ii}^A$), i.e., as weighted averages of all the poles of $G_{ii}^R$ ($G_{ii}^A$) within the basis that diagonalizes the Λ^R (Λ^A) matrix. Indeed, inserting a complete set of eigenstates of the (N − 1)-electron system in Eq. 7, the commutator can be rewritten as

$\begin{array}{ccc}\hfill {ϵ}_i^R& =\hfill & -\sum _k\frac{\langle {\mathrm{\Psi }}_0^N|{\widehat{O}}_i^{†}|{\mathrm{\Psi }}_k^{N-1}\rangle \langle {\mathrm{\Psi }}_k^{N-1}|[\widehat{H},{\widehat{O}}_i]|{\mathrm{\Psi }}_0^N\rangle }{\langle {\mathrm{\Psi }}_0^N|{\widehat{O}}_i^{†}{\widehat{O}}_i|{\mathrm{\Psi }}_0^N\rangle }\hfill \\ \hfill & =\hfill & -\sum _k\frac{\langle {\mathrm{\Psi }}_0^N|{\widehat{O}}_i^{†}|{\mathrm{\Psi }}_k^{N-1}\rangle \langle {\mathrm{\Psi }}_k^{N-1}|{\widehat{O}}_i|{\mathrm{\Psi }}_0^N\rangle }{\langle {\mathrm{\Psi }}_0^N|{\widehat{O}}_i^{†}{\widehat{O}}_i|{\mathrm{\Psi }}_0^N\rangle }(E_k^{N-1}-E_0^N)\hfill \\ \hfill & =\hfill & \frac{{\sum }_k{B}_{ii}^{k,R}{ϵ}_k^R}{{\sum }_k{B}_{ii}^{k,R}},\hfill \end{array}(16)$

which is a weighted average of the poles of $G_{ii}^R$. Inserting a complete set of eigenstates of the (N + 1)-electron system in Eq. 10 one can show in a similar way that the ${ϵ}_i^A={\sum }_k{B}_{ii}^{k,A}{ϵ}_k^A/{\sum }_k{B}_{ii}^{k,A}$ within the EKT basis. This means that if there are not satellites in the EKT basis set, then the EKT removal/addition energies are exact, provided that one uses the exact first moment.

We notice that very recently Lee et al. (Lee et al., 2021), have also proposed an expression for the spectral function from the EKT.

3.1.2 Insights Into the MEET

Several choices for an optimal basis set for the MEET expressions are now possible. In previous works we considered the basis of natural orbitals as optimal basis set for the MEET based on the following findings (Di Sabatino, 2016): i) the MEET (in this basis of natural orbitals) gives the exact spectral function at all level of approximations for the symmetric Hubbard dimer using exact density matrices; ii) the MEET in its lowest level of approximation in terms of one- and two-body density matrices gives good results for the spectral function of the (symmetric) Hubbard model with more sites using approximate density matrices. Moreover, for these (symmetric) model systems the Λ^R/A matrices of the EKT are diagonal in the basis of natural orbitals, therefore there is not another better option. For the asymmetric Hubbard dimer instead, and in general for realistic systems, the basis of natural orbitals does not diagonalize the Λ^R/A matrices, therefore the set which diagonalizes these matrices can be a better option for the MEET. We notice that this choice of the optimal basis set can be generally applied to other methods which express G as a continued fraction, such as the Lanczos method (Balzer et al., 2011), in order to have more accurate results at a given order of truncation of the series. For example in Figure 1 we report the spectral function of the Hubbard dimer governed by the Hamiltonian in Eq. 13 for two different values of |U₁ − U₂|. The (D)EKT results are obtained using exact density matrices. The results show that the basis which diagonalizes the Λ^R/A matrices is a much better choice than the basis of natural orbitals the more the difference |U₁ − U₂| is large. We also observe that the removal part is less affected by the diagonal approximation than the addition part, and we observe this trend also in more complex systems. The diagonal approximation has been addressed in literature also for realistic systems. (Piris et al., 2013; Kent et al., 1998) In particular in bulk silicon QMC results show that the DEKT slightly overestimates the EKT band gaps. Below we will address this system in more details.

FIGURE 1

FIGURE 1. Spectral function of the Hubbard dimer with two different on-site interactions U₁ and U₂ for site 1 and site 2, respectively: exact solution (black curves) vs. EKT (black curves, EKT is exact in this case, see text) and DEKT/MEET in the basis of natural orbitals (red curves).

3.2 Realistic Systems: The Example of Bulk Si

As an example of realistic systems we use bulk silicon, for which results using the EKT within QMC are reported in Ref. (Kent et al., 1998). The diagonal approximation to the EKT within QMC works very well for the valence states and slightly less well for the conduction states, with a band gap at the Γ point of 4.4 eV vs. 3.8 eV from the full EKT compared to 3.4 eV in experiment. (Kent et al., 1998) We note that this discrepency is largely due to the energy of the conduction band at Γ. However, bulk silicon is a relatively weakly correlated system, with hence a predominant quasiparticle-like spectral function for which the EKT is a good approximation. Larger overestimation of the band gap can be expected for strongly correlated systems. This can be understood from the interpretation of the EKT energies as first moments of the one-body Green’s function. However, an important point to stress is that even for these systems, which are a challenge for state-of-the-art ab initio methods, such as GW, the EKT would open a gap, in accordance with experiment.

As pointed out in Ref. (Kent et al., 1998) the choice of the trial wave function for QMC calculations is of critical importance. Indeed, as a result of the fixed-node approximation, QMC calculations of the matrix elements of the density matrix and operators V^R,A (see EKT Equations 8, 11) are expected to critically depend on the nodal structure of the trial wave function employed. For the weakly correlated bulk silicon the accurate QMC value of 3.8 eV reported above has been obtained by (Kent et al., 1998) using a standard Slater-Jastrow trial wave function whose nodes are those of a single determinant consisting of LDA orbitals. For more strongly correlated systems the wave function acquires a significant multi-determinant character and getting physically meaningful nodes becomes much more difficult. It is thus useful to use the EKT within alternative approaches.

In our previous works (Di Sabatino et al., 2016, 2019, 2021) we used reduced-density matrix functional theory (RDMFT) (Gilbert, 1975) to find approximations to the one- and two-body density matrices which are needed in the MEET equations. More specifically the two-body density matrix in the MEET equations is approximated using the Power functional, which is given by ${\mathrm{\Gamma }}_{ijkl}^{(2)}=n_i{n}_j{\delta }_{il}{\delta }_{jk}-n_i^{\alpha }{n}_j^{\alpha }{\delta }_{ik}{\delta }_{jl}$ (α = 0.65). (Sharma et al., 2013) The optimal natural orbitals {ϕ_i} and occupation numbers {n_i} are obtained by minimizing the total energy which is expressed in terms of γ and Γ⁽²⁾, with Γ⁽²⁾ as functional of γ. The Power functional is used also to approximate Γ⁽²⁾[γ] in the energy functional. In this work we use the same protocol for the EKT equations. We implemented the EKT equations in a modified version of the full-potential linearized augmented plane wave (FP-LAPW) code ELK (Elk, 2004), with practical details of the calculations following the scheme described in Ref. (Sharma et al., 2008). For bulk Si we used a lattice constant of 5.43 Å  and a Γ-centered 8 × 8 × 8 k-point sampling of the Brillouin zone. In Figure 2 we report the DEKT spectral function of bulk silicon: the direct band gap at Γ is 12.9 eV, while the fundamental band gap is 8.18 eV, which is larger than the experimental one of 1.12 eV (Sze, 1969). We also observe a spurious peak in the band gap due to the fact that the Power functional produces occupation numbers which strongly deviate from 1 and 0 (as one would expect for this weakly correlated system) close to the Fermi energy (see bottom panel of Figure 2). This is in contrast with the QMC results. Note that we observe similar deviations from 1 and 0 also for other weakly correlated systems, such as diamond, which points to a problem of the Power functional for the description of occupation numbers of weakly correlated systems. Moreover, the full EKT does not show any improvement over the DEKT, as one can see from Figure 3, in which the EKT and DEKT energies are reported: the fundamental band gap is reduced by only 0.06 eV. This is again in contrast with the QMC results in which, although small, there is a significant difference. We attribute this different trend to the use of the Power functional, which contracts the four-point 2-RDM to two points only, and hence probably mitigating the impact of the diagonalization of the Λ^R/A matrices. These results on bulk Si indicate that, although the EKT/DEKT are expected to overestimate the band gap (even using very accurate density matrices, as for example shown in the case of the Hubbard model (Di Sabatino et al., 2016)), this overestimation can be much amplified by using approximations such as the Power functional. More advanced approximations to Γ⁽²⁾ are hence needed, which give, in particular, more accurate natural occupation numbers. We notice that varying α would change the band gap width. In particular α = 1 would give the HF band gap, which still overestimates the experimental one, whereas decreasing α would increase the overestimation of the band gap.

FIGURE 2

FIGURE 2. Spectral function of bulk Si within DEKT (violet solid line). The experimental photoemission spectrum (small triangles) is taken from Ref. (Chelikowsky et al., 1989). The color map illustrates the occupation numbers n_i that play a role into the spectrum for the reported energy range.

FIGURE 3

FIGURE 3. Removal/addition energies of bulk Si: EKT (filled circles) vs DEKT (empty circles) for two k-points, k = (0, 0, 0) and k = (3/4, 1/4, 0).

4 Conclusions and Perspectives

We linked our recently derived Many-Body Effective Energy Theory (MEET) for the calculation of photoemission spectra to the Extended Koopman’s Theorem (EKT). Within the lowest level of approximation in terms of one- and two-body density matrices, the MEET equations correspond to the so-called diagonal approximation to the EKT (DEKT) equations. This allowed us to readily extend the EKT to the calculation of an approximate spectral function as well as to give an alternative interpretation of the EKT in terms of moments of the one-body Green’s function. Using the test case of the Hubbard dimer with two different on-site interactions U₁ and U₂ for site 1 and site 2 we showed the effect of the basis set on the MEET (removal and addition) energies: in particular HOMO-LUMO gap in the basis sets which solve the EKT secular equations (one basis set for the valence part and one for the conduction part) is smaller than the HOMO-LUMO gap obtained using the natural orbital basis set. These results are in line with the EKT results reported in literature for bulk Si using QMC. We have implemented the EKT within reduced-density matrix functional theory (RDMFT), which offers a convenient computationally affordable framework to treat extended systems. However one has to rely on approximate one- and two-body density matrices. We showed that using the currently available approximations the DEKT band gap of Si largely deviates from the DEKT value obtained using QMC (12.9 eV vs 4.4 eV at the Γ point) and, moreover, there is no effect of the basis set (EKT vs DEKT) on the DEKT energies, contrary to what is observed within QMC, where, although small, there is a significant difference. These results on bulk Si indicate that, although the EKT/DEKT are expected to overestimate the band gap (even using very accurate density matrices), this overestimation can be much amplified by commonly used approximations in RDMFT. This also explains the huge overestimation of the band gap obtained by the MEET within RDMFT in strongly correlated systems such as paramagnetic NiO. We are currently working on improving approximations to correlation in RDMFT by introducing some form of screening (for example the screening due to electron-hole excitations as in GW), which is of particular importance in solids.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author Contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

This study has been supported through the EUR grant NanoX ANR-17-EURE-0009 in the framework of the “Programme des Investissements d’Avenir” and by ANR (project ANR-18-CE30-0025 and ANR-19-CE30-0011).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

¹Since in the basis of natural orbitals the S^R (S^A) matrix is a diagonal matrix with the natural occupation numbers n_i (1-n_i) as elements, the invertibility of this matrix is strictly related to the non-existence of so-called pinned states, i.e. states with occupation numbers equal to 1 or 0. This is an important question that has several consequences (e.g., Giesbertz and van Leeuwen, 2013; Baldsiefen et al., 2015). Here we assume that S^R (S^A) is invertible in a restricted space (of natural orbitals) in which the corresponding KS orbitals are occupied (unoccupied). This is a reasonable assumption.

²Note that the standard EKT equations in RDMFT (Morrell et al., 1975) present a prefactor $1/\sqrt{n_i{n}_j}$ ($1/\sqrt{(1-n_i)(1-n_j)}$) instead of 1/n_i (1/(1 − n_i)) in Eq. 9 (Eq. 12); these two choices yield the same eigenvalues for the matrix Λ^R (Λ^A).

References

Baldsiefen, T., Cangi, A., and Gross, E. K. U. (2015). Reduced-density-matrix-functional Theory at Finite Temperature: Theoretical Foundations. Phys. Rev. A. 92, 052514. doi:10.1103/PhysRevA.92.052514

CrossRef Full Text | Google Scholar

Balzer, M., Gdaniec, N., and Potthoff, M. (2011). Krylov-space Approach to the Equilibrium and Nonequilibrium Single-Particle Green's Function. J. Phys. Condens. Matter 24, 035603. doi:10.1088/0953-8984/24/3/035603

CrossRef Full Text | Google Scholar

Berger, J. A., Romaniello, P., Tandetzky, F., Mendoza, B. S., Brouder, C., and Reining, L. (2014). Solution to the many-body Problem in One point. New J. Phys. 16, 113025. doi:10.1088/1367-2630/16/11/113025

CrossRef Full Text | Google Scholar

Bozkaya, U. (2013). The Extended Koopmans' Theorem for Orbital-Optimized Methods: Accurate Computation of Ionization Potentials. J. Chem. Phys. 139, 154105. doi:10.1063/1.4825041

CrossRef Full Text | Google Scholar

Bozkaya, U., and Ünal, A. (2018). State-of-the-Art Computations of Vertical Ionization Potentials with the Extended Koopmans' Theorem Integrated with the CCSD(T) Method. J. Phys. Chem. A. 122, 4375–4380. doi:10.1021/acs.jpca.8b01851

PubMed Abstract | CrossRef Full Text | Google Scholar

Chelikowsky, J. R., Wagener, T. J., Weaver, J. H., and Jin, A. (1989). Valence-and Conduction-Band Densities of States for Tetrahedral Semiconductors: Theory and Experiment. Phys. Rev. 40, 9644–9651. doi:10.1103/PhysRevB.40.9644

PubMed Abstract | CrossRef Full Text | Google Scholar

Cioslowski, J., Piskorz, P., and Liu, G. (1997). Ionization Potentials and Electron Affinities from the Extended Koopmans' Theorem Applied to Energy-Derivative Density Matrices: The EKTMPn and EKTQCISD Methods. J. Chem. Phys. 107, 6804–6811. doi:10.1063/1.474921

CrossRef Full Text | Google Scholar

Dahlen, N. E., and van Leeuwen, R. (2005). Self-consistent Solution of the Dyson Equation for Atoms and Molecules within a Conserving Approximation. J. Chem. Phys. 122, 164102. doi:10.1063/1.1884965

CrossRef Full Text | Google Scholar

Di Sabatino, S., Berger, J. A., Reining, L., and Romaniello, P. (2016). Photoemission Spectra from Reduced Density Matrices: The Band gap in Strongly Correlated Systems. Phys. Rev. B 94, 155141. doi:10.1103/PhysRevB.94.155141

CrossRef Full Text | Google Scholar

Di Sabatino, S., Berger, J. A., Reining, L., and Romaniello, P. (2015). Reduced Density-Matrix Functional Theory: Correlation and Spectroscopy. J. Chem. Phys. 143, 024108. doi:10.1063/1.4926327

CrossRef Full Text | Google Scholar

Di Sabatino, S., Berger, J. A., and Romaniello, P. (2019). Many-body Effective Energy Theory: Photoemission at strong Correlation. J. Chem. Theor. Comput. 15, 5080–5086. doi:10.1021/acs.jctc.9b00427

CrossRef Full Text | Google Scholar

Di Sabatino, S., Koskelo, J., Berger, J. A., and Romaniello, P. (2021). Photoemission Spectrum in Paramagnetic Feo under Pressure: Towards an Ab Initio Description. Phys. Rev. Res. 3, 013172. doi:10.1103/PhysRevResearch.3.013172

CrossRef Full Text | Google Scholar

Di Sabatino, S. (2016). Private Communication.

Google Scholar

Elk (2004). Elk. Available at: http://elk.sourceforge.net.

Google Scholar

Giesbertz, K. J. H., and van Leeuwen, R. (2013). Natural Occupation Numbers: When Do They Vanish? J. Chem. Phys. 139, 104109. doi:10.1063/1.4820419

CrossRef Full Text | Google Scholar

Gilbert, T. L. (1975). Hohenberg-kohn Theorem for Nonlocal External Potentials. Phys. Rev. B 12, 2111–2120. doi:10.1103/PhysRevB.12.2111

CrossRef Full Text | Google Scholar

Guzzo, M., Lani, G., Sottile, F., Romaniello, P., Gatti, M., Kas, J. J., et al. (2011). Valence Electron Photoemission Spectrum of Semiconductors:Ab InitioDescription of Multiple Satellites. Phys. Rev. Lett. 107, 166401. doi:10.1103/physrevlett.107.166401

PubMed Abstract | CrossRef Full Text | Google Scholar

Hedin, L. (1965). New Method for Calculating the One-Particle Green’s Function with Application to the Electron-Gas Problem. Phys. Rev. 139, A796–A823. doi:10.1103/physrev.139.a796

CrossRef Full Text | Google Scholar

Katriel, J., and Davidson, E. R. (1980). Asymptotic Behavior of Atomic and Molecular Wave Functions. Proc. Natl. Acad. Sci. 77, 4403–4406. doi:10.1073/pnas.77.8.4403

PubMed Abstract | CrossRef Full Text | Google Scholar

Kent, P. R. C., Hood, R. Q., Towler, M. D., Needs, R. J., and Rajagopal, G. (1998). Quantum Monte Carlo Calculations of the One-Body Density Matrix and Excitation Energies of Silicon. Phys. Rev. B 57, 15293–15302. doi:10.1103/PhysRevB.57.15293

CrossRef Full Text | Google Scholar

Kuneš, J., Anisimov, V. I., Skornyakov, S. L., Lukoyanov, A. V., and Vollhardt, D. (2007). Nio: Correlated Band Structure of a Charge-Transfer Insulator. Phys. Rev. Lett. 99, 156404. doi:10.1103/PhysRevLett.99.156404

PubMed Abstract | CrossRef Full Text | Google Scholar

Lani, G., Romaniello, P., and Reining, L. (2012). Approximations for many-body Green's Functions: Insights from the Fundamental Equations. New J. Phys. 14, 013056. doi:10.1088/1367-2630/14/1/013056

CrossRef Full Text | Google Scholar

Lee, J., Malone, F. D., Morales, M. A., and Reichman, D. R. (2021). Spectral Functions from Auxiliary-Field Quantum Monte Carlo without Analytic Continuation: The Extended Koopmans' Theorem Approach. J. Chem. Theor. Comput. 17, 3372–3387. PMID: 33983735. doi:10.1021/acs.jctc.1c00100

CrossRef Full Text | Google Scholar

Leiva, P., and Piris, M. (2005). Assessment of a New Approach for the Two-Electron Cumulant in Natural-Orbital-Functional Theory. J. Chem. Phys. 123, 214102. doi:10.1063/1.2135289

CrossRef Full Text | Google Scholar

Lischner, J., Vigil-Fowler, D., and Louie, S. G. (2013). Physical Origin of Satellites in Photoemission of Doped Graphene: AnAb InitioGWPlus Cumulant Study. Phys. Rev. Lett. 110, 146801. doi:10.1103/PhysRevLett.110.146801

PubMed Abstract | CrossRef Full Text | Google Scholar

Morrell, M. M., Parr, R. G., and Levy, M. (1975). Calculation of Ionization Potentials from Density Matrices and Natural Functions, and the Long-Range Behavior of Natural Orbitals and Electron Density. J. Chem. Phys. 62, 549–554. doi:10.1063/1.430509

CrossRef Full Text | Google Scholar

Pavlyukh, Y. (2019). The Ubiquitous Extended Koopmans' Theorem. Phys. Status Solidi B 256, 1800591. doi:10.1002/pssb.201800591

CrossRef Full Text | Google Scholar

Pernal, K., and Cioslowski, J. (2005). Ionization Potentials from the Extended Koopmans' Theorem Applied to Density Matrix Functional Theory. Chem. Phys. Lett. 412, 71–75. doi:10.1016/j.cplett.2005.06.103

CrossRef Full Text | Google Scholar

Piris, M., Matxain, J. M., Lopez, X., and Ugalde, J. M. (2012). The Extended Koopmans’ Theorem: Vertical Ionization Potentials from Natural Orbital Functional Theory. J. Chem. Phys. 136, 174116. doi:10.1063/1.4709769

CrossRef Full Text | Google Scholar

Piris, M., Matxain, J. M., Lopez, X., and Ugalde, J. M. (2013). The One-Electron Picture in the Piris Natural Orbital Functional 5 (Pnof5). Theor. Chem. Acc. 132, 1298. doi:10.1007/s00214-012-1298-4

CrossRef Full Text | Google Scholar

Romaniello, P., Bechstedt, F., and Reining, L. (2012). Beyond theGWapproximation: Combining Correlation Channels. Phys. Rev. B 85, 155131. doi:10.1103/physrevb.85.155131

CrossRef Full Text | Google Scholar

Romaniello, P., Guyot, S., and Reining, L. (2009). The Self-Energy beyond Gw: Local and Nonlocal Vertex Corrections. J. Chem. Phys. 131, 154111. doi:10.1063/1.3249965

CrossRef Full Text | Google Scholar

Sharma, S., Dewhurst, J. K., Lathiotakis, N. N., and Gross, E. K. U. (2008). Reduced Density Matrix Functional for many-electron Systems. Phys. Rev. B 78, 201103. doi:10.1103/PhysRevB.78.201103

CrossRef Full Text | Google Scholar

Sharma, S., Dewhurst, J. K., Shallcross, S., and Gross, E. K. U. (2013). Spectral Density and Metal-Insulator Phase Transition in mott Insulators within Reduced Density Matrix Functional Theory. Phys. Rev. Lett. 110, 116403. doi:10.1103/physrevlett.110.116403

PubMed Abstract | CrossRef Full Text | Google Scholar

Shishkin, M., Marsman, M., and Kresse, G. (2007). Accurate Quasiparticle Spectra from Self-ConsistentGWCalculations with Vertex Corrections. Phys. Rev. Lett. 99, 246403. doi:10.1103/PhysRevLett.99.246403

PubMed Abstract | CrossRef Full Text | Google Scholar

Smith, D. W., and Day, O. W. (1975). Extension of Koopmans' Theorem. I. Derivation. J. Chem. Phys. 62, 113–114. doi:10.1063/1.430253

CrossRef Full Text | Google Scholar

Springer, M., Aryasetiawan, F., and Karlsson, K. (1998). First-PrinciplesT-Matrix Theory with Application to the 6 eV Satellite in Ni. Phys. Rev. Lett. 80, 2389–2392. doi:10.1103/PhysRevLett.80.2389

CrossRef Full Text | Google Scholar

Stan, A., Dahlen, N. E., and Leeuwen, R. v. (2006). Fully Self-Consistent GW Calculations for Atoms and Molecules. Europhys. Lett. 76, 298–304. doi:10.1209/epl/i2006-10266-6

CrossRef Full Text | Google Scholar

Stan, A., Dahlen, N. E., and van Leeuwen, R. (2009). Levels of Self-Consistency in the Gw Approximation. J. Chem. Phys. 130, 114105. doi:10.1063/1.3089567

CrossRef Full Text | Google Scholar

Stefanucci, G., Pavlyukh, Y., Uimonen, A.-M., and van Leeuwen, R. (2014). Diagrammatic Expansion for Positive Spectral Functions beyondGW: Application to Vertex Corrections in the Electron Gas. Phys. Rev. B 90, 115134. doi:10.1103/PhysRevB.90.115134

CrossRef Full Text | Google Scholar

Sundholm, D., and Olsen, J. (1993). The Exactness of the Extended Koopmans' Theorem: A Numerical Study. J. Chem. Phys. 98, 3999–4002. doi:10.1063/1.464028

CrossRef Full Text | Google Scholar

Sze, S. M. (1969). The Physics of Semiconductor Devices. New York: Wiley.

Google Scholar

Zheng, H. (2016). First Principles Quantum Monte Carlo Study of Correlated Electronic Systems. University of Illinois at Urbana-Champaign. Ph.D. thesis.

Google Scholar

Zhukov, V. P., Chulkov, E. V., and Echenique, P. M. (2004). Lifetimes of Excited Electrons in Fe and Ni: First-Principles GW and theT-Matrix Theory. Phys. Rev. Lett. 93, 096401. doi:10.1103/PhysRevLett.93.096401

PubMed Abstract | CrossRef Full Text | Google Scholar