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 δ i R / A ( ω ) as (Di Sabatino et al., 2016): G i i R ( ω ) = γ i i ω − δ i R ( ω ) − i η , (1) G i i A ( ω ) = 1 − γ i i ω − δ i A ( ω ) + i η , (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 i i ( ω ) = γ i i δ ( ω − δ i R ( ω ) ) + ( 1 − γ i i ) δ ( ω − δ i A ( ω ) ) , (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 δ i R / A ( ω ) 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). δ i R , ( 1 ) = h i i + 1 n i ∑ k l m V i m k l Γ k l m i ( 2 ) (4) δ i A , ( 1 ) = h i i + 1 ( 1 − n i ) ∑ k V i k i k − V i k k i n k − 1 ( 1 − n i ) ∑ k l m V i m k l Γ k l m i ( 2 ) , (5) where Γ k l m i ( 2 ) = ⟨ Ψ 0 N | a ̂ i † a ̂ m † a ̂ l a ̂ k | Ψ 0 N ⟩ are the matrix elements of the two-body reduced density matrix, with Ψ 0 N the ground-state wavefunction of the N-electron system. As discussed in Ref. (Di Sabatino et al., 2016) the various approximations δ i R / A , ( n ) ( ω ) are related to the n-th moments μ n , i R / A = ∑ k B i i k , R / A ( ϵ k R / A ) n ∑ k B i i k , R / A of the G i i R / A ( ω ) . 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 B i i k , R = 〈 Ψ 0 N | c ̂ i † | Ψ k N − 1 〉 〈 Ψ k N − 1 | c ̂ i | Ψ 0 N 〉 B i i k , A = 〈 Ψ 0 N | c ̂ i | Ψ k N + 1 〉 〈 Ψ k N + 1 | c ̂ i † | Ψ 0 N 〉 , with E 0 N and Ψ 0 N the ground-state energy and wave function of the N-electron system and E k N ± 1 and Ψ k N ± 1 the kth state energy and wave function of the (N ± 1)-electron system. This allows for a more compact expression of G i i R / A ( ω ) as a continued fraction of moments G i i R = n i ω − μ 1 , i R ω − μ 1 , i R … ω − μ 2 , i R μ 1 , i R … , (6) (and similarly for G i i 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 i i R / A has only one pole which is a weighted average of all the poles of G i i 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.

Within the EKT one starts from the following approximation for the removal energy ϵ i R (Kent et al., 1998) ϵ i R = − 〈 Ψ 0 N | O ̂ i † [ H ̂ , O ̂ i ] | Ψ 0 N 〉 〈 Ψ 0 N | O ̂ i † O ̂ i | Ψ 0 N 〉 (7) with Ψ 0 N the ground-state many-body wave function of the N-electron system, and O ̂ i = ∑ k C k i R a ̂ k , O ̂ i † = ∑ k C k i R * a ̂ k † , with { C k i R } a set of coefficients to be determined. The stationary condition (with respect to the coefficients C k i R ) for ϵ i R leads to the secular equation ( V R − ϵ i R S R ) C i R = 0 , (8) with V i j R = − ⟨ Ψ 0 N | a ̂ j † [ H ̂ , a ̂ i ] | Ψ 0 N ⟩ and S ^R the one-body density matrix S i j R = γ i j = ⟨ Ψ 0 N | a ̂ j † a ̂ i | Ψ 0 N ⟩ . If one defines the matrix Λ R = [ S R ] − 1 V R in the basis of natural orbitals, with S i j R = n i δ i j and works out the commutator in V i j R , one arrives at Λ i j R = 1 n i n i h j i + ∑ k l m V j m k l Γ k l m i ( 2 ) . (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 = 〈 Ψ 0 N | [ H ̂ , O ̂ i ] O ̂ i † | Ψ 0 N 〉 〈 Ψ 0 N | O ̂ i O ̂ i † | Ψ 0 N 〉 (10) and in a similar way as for ϵ i R we arrive at the eigenvalue equation ( V A − ϵ i A S A ) C i A = 0 , (11) with V i j A = ⟨ Ψ 0 N | a ̂ i [ H ̂ , a ̂ j † ] | Ψ 0 N ⟩ and S ^A related to the one-body density matrix as S i j A = 1 − γ i j . Similarly to the removal energy problem, using the basis of natural orbitals, one can work out the commutator in V i j A and reformulate the problem in terms of the matrix Λ A = [ S A ] − 1 V A ¹ , which reads Λ i j A = 1 ( 1 − n i ) × ( 1 − n i ) h j i + ∑ k V j k i k − V j k k i n k − ∑ k l m V j m k l Γ k l m i ( 2 ) . (12) Again, the diagonal elements of Λ ^A are the MEET addition energies within the approximation given in Eq. 5. ²

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 ∑ i , j = 1,2 i ≠ j ∑ σ a ̂ i σ † a ̂ j σ + U 1 n ̂ 1 ↑ n ̂ 1 ↓ + U 2 n ̂ 2 ↑ n ̂ 2 ↓ , (13) where i, j run over the sites, n ̂ i σ = a ̂ i σ † a ̂ i σ , 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.

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 Γ i j k l ( 2 ) = n i n j δ i l δ j k − n i α n j α δ i k δ j l (α = 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.