Multiple impurities and combined local density approximations in Site-Occupation Embedding Theory

Multiple impurities and combined local density approximations in Site-Occupation Embedding Theory

Bruno Senjean Laboratoire de Chimie Quantique, Institut de Chimie, CNRS/Université de Strasbourg, 4 rue Blaise Pascal, 67000 Strasbourg, France    Naoki Nakatani Department of Chemistry, Graduate School of Science and Engineering, Tokyo Metropolitan University, 1-1 Minami-Osawa, Hachioji, Tokyo 192-0397, Japan Institute for Catalysis, Hokkaido University, N21W10 Kita-ku, Sapporo, Hokkaido 001-0021, Japan    Masahisa Tsuchiizu Department of Physics, Nara Women’s University, Nara 630-8506, Japan    Emmanuel Fromager Laboratoire de Chimie Quantique, Institut de Chimie, CNRS/Université de Strasbourg, 4 rue Blaise Pascal, 67000 Strasbourg, France

Site-occupation embedding theory (SOET) is an in-principle-exact multi-determinantal extension of density-functional theory for model Hamiltonians. Various extensions of recent developments in SOET [Senjean et al., Phys. Rev. B 97, 235105 (2018)] are explored in this work. An important step forward is the generalization of the theory to multiple impurity sites. We also propose a new single-impurity density-functional approximation (DFA) where the density-functional impurity correlation energy of the two-level (2L) Hubbard system is combined with the Bethe ansatz local density approximation (BALDA) to the full correlation energy of the (infinite) Hubbard model. In order to test the new DFAs, the impurity-interacting wavefunction has been computed self-consistently with the density matrix renormalization group method (DMRG). Double occupation and per-site energy expressions have been derived and implemented in the one-dimensional case. A detailed analysis of the results is presented, with a particular focus on the errors induced either by the energy functionals solely or by the self-consistently converged densities. Among all the DFAs (including those previously proposed), the combined 2L-BALDA is the one that performs the best in all correlation and density regimes. Finally, extensions in new directions, like a partition-DFT-type reformulation of SOET, a projection-based SOET approach, or the combination of SOET with Green functions, are briefly discussed as a perspective.

thanks: Corresponding author

I Introduction

The accurate and low-cost description of strongly correlated materials remains one of the most challenging task in electronic structure theory As highly accurate wavefunction-based methods are too expensive to be applied to the whole system of interest, simplified and faster solutions have to be considered. Such simplications should ideally not alter the description of strong correlation effects. This is where the challenge stands. Based on the cogent argument that strong electron correlation is essentially local Pulay (1983); Saebo and Pulay (1993); Hampel and Werner (1996) and that the region of interest is one part of a much larger (extended) system, embedding approaches are mainly used in practice Sun and Chan (2016). The basic idea is to map the fully-interacting problem onto a so-called impurity-interacting one. In the Hubbard model, the impurity corresponds to an atomic site. Among such embedding techniques are the well-established dynamical mean-field theory (DMFT) Georges and Kotliar (1992); Georges et al. (1996); Kotliar and Vollhardt (2004); Held (2007); Zgid and Chan (2011), its cluster  Hettler et al. (1998, 2000); Lichtenstein and Katsnelson (2000); Kotliar et al. (2001); Maier et al. (2005) and diagrammatic Rohringer et al. (2018) extensions, as well as combinations of DMFT with either density-functional theory (DFT) [the so-called DMFT+DFT approach Kotliar et al. (2006)] or the Green-function-based GW method  Sun and Kotliar (2002); Biermann et al. (2003); Karlsson (2005); Boehnke et al. (2016); Werner and Casula (2016); Nilsson et al. (2017). Such combinations aim at incorporating non-local correlation effects in DMFT. More recently, the self-energy embedding theory (SEET) Kananenka et al. (2015); Lan et al. (2015); Zgid and Gull (2017); Lan et al. (2017), which can be applied to both model and ab initio Hamiltonians, has been developed. Let us stress that all the aforementioned embedding techniques use the frequency-dependent one-particle Green function as basic variable.

Alternative frequency-independent approaches like the density-matrix embedding theory (DMET) Knizia and Chan (2012, 2013); Bulik et al. (2014); Zheng and Chan (2016); Wouters et al. (2016, 2017); Rubin (2016) have emerged in recent years. By construction, standard approximate DMET does not describe correlation effects in the environment, thus requiring the treatment of more than one impurity site in order to obtain reasonably accurate results Knizia and Chan (2012). More correlation can be incorporated into DMET by using an antisymmetrized geminal power wavefunction Tsuchimochi et al. (2015), or, alternatively, by improving the description of the boundary between the fragment and the bath Welborn et al. (2016). Note that Ayral et al. Ayral et al. (2017) succeeded recently in establishing formal connections between DMET, DMFT and the rotationally invariant slave bosons (RISB) theory Frésard and Wölfle (1992); Lechermann et al. (2007).

This paper deals with another frequency-independent approach, namely site-occupation embedding theory (SOET) Fromager (2015); Senjean et al. (2017, 2018). While, in conventional Kohn–Sham (KS) DFT, the fully-interacting problem is mapped onto a non-interacting one, an auxiliary impurity-interacting system is used in SOET for extracting the density (i.e. the site occupations in this context) and, through an appropriate density functional for the environment (referred to as bath), the total energy. From a quantum chemical point of view, SOET is nothing but a multi-determinantal extension of KS-DFT for model Hamiltonians Chayes et al. (1985); Gunnarsson and Schönhammer (1986); Capelle and Campo Jr. (2013). In a recent paper Senjean et al. (2018), the authors explained how exact expressions for the double occupation and the per-site energy of the uniform Hubbard model can be extracted from SOET. They also proposed various local density-functional approximations for the bath. The latter work suffered from two main weaknesses. First of all, the complete self-consistent formulation of the theory was done only for a single impurity site, thus preventing a gradual transition from KS-DFT (no impurity sites) to pure wavefunction theory (no bath sites). Moreover, none of the proposed DFAs gave satisfactory results in all correlation and density regimes.

We explain in this work how these limitations can be overcome. The paper is organized as follows. First, an in-principle-exact generalization of SOET to multiple impurity sites is derived in Sec. II.1. The resulting expressions for the double occupation and the per-site energy in the uniform case are given in Sec. II.2. Existing and newly proposed DFAs are then discussed in detail in Sec. III. Following the computational details in Sec. IV, results obtained at half-filling (Sec. V.1) and away from half-filling (Secs. V.2 and V.3) are presented and analyzed. Exact properties of the impurity correlation potential are discussed in Sec. V.4. Conclusions and perspectives are finally given in Sec. VI.

Ii Theory

ii.1 Site-occupation embedding theory with multiple impurities

Let us consider the (not necessarily uniform) -site Hubbard Hamiltonian with external potential ,


The hopping operator, which is the analog for model Hamiltonians of the kinetic energy operator, reads as follows in second quantization,


where is the so-called hopping parameter and means that the atomic sites and are nearest neighbors. The on-site two-electron repulsion operator with strength and the local external potential operator (which is the analog for model Hamiltonians of the nuclear potential) are expressed in terms of the spin-density and density operators as follows,





The exact ground-state energy of can be obtained variationally as follows, in complete analogy with conventional DFT Hohenberg and Kohn (1964),


where is a trial collection of site occupations (simply called density in the following) and . Within the Levy–Lieb (LL) constrained-search formalism Levy (1979), the Hohenberg–Kohn functional can be rewritten as follows in this context,


where the minimization is restricted to wavefunctions with density . As shown in previous works Fromager (2015); Senjean et al. (2017, 2018), the exact minimizing density in Eq. (5) can be obtained from a fictitious partially-interacting system consisting of interacting impurity sites surrounded by non-interacting ones (the so-called bath sites), thus leading to an in-principle-exact SOET. While our recent developments focused on the single-impurity version of SOET, we propose in the following a general formulation of the theory with an arbitrary number of impurity sites. Such a formulation was briefly mentioned in Ref. Fromager (2015) for the purpose of deriving an adiabatic connection formula for the correlation energy of the bath.

Let us introduce the analog for impurity sites of the LL functional,


where . Note that, for convenience, the impurity sites have been labelled as . If we now introduce the complementary Hartree-exchange-correlation (Hxc) functional for the bath,


the ground-state energy expression in Eq. (5) can be rewritten as follows,


or, equivalently,


where , thus leading to the final variational expression

The minimizing -impurity-interacting wavefunction in Eq. (II.1) reproduces the exact density profile of the fully-interacting system described by the Hubbard Hamiltonian in Eq. (1). From the stationarity in of the energy, we obtain the following self-consistent equation,




plays the role of an embedding potential for the impurities. In the particular case of a uniform half-filled density profile, the embedding potential equals zero in the bath and on the impurity sites. This key result, which appears when applying the hole-particle symmetry transformation to the impurity-interacting LL functional, is proved in Appendix A, thus providing a generalization of Appendix C in Ref. Senjean et al. (2018). Note that the KS and Schrödinger equations are recovered from Eq. (12) when and (i.e. the total number of sites), respectively. In SOET, is in the range , thus leading to a hybrid formalism where a many-body correlated wavefunction is embedded into a DFT potential. In practice, Eq. (12) can be solved, for example, by applying an exact diagonalization procedure Senjean et al. (2017) (which corresponds to a full configuration interaction) for small rings, or by using the more advanced density matrix renormalization group (DMRG) method which allows for the description of larger systems Senjean et al. (2018).

Let us now return to the expression in Eq. (8) of the complementary Hxc energy for the bath. By using the KS decompositions of the fully-interacting and -impurity-interacting LL functionals,




respectively, where the non-interacting kinetic energy functional reads as follows in this context,


we obtain


If we now separate the Hxc energies into Hx (i.e. mean-field) and correlation contributions,




we obtain the final expression




While local density approximations (LDA) based, for example, on the Bethe ansatz (BALDA) are available for in the literature Lima et al. (2002, 2003); Capelle et al. (2003), no DFA has been developed so far for modeling the correlation energy of multiple impurity sites. Existing approximations for a single impurity are reviewed in Sec. III.2. Newly proposed DFAs will be introduced in Secs. III.3 and III.4.

ii.2 Exact double occupation and per-site energy expressions in the uniform case

In this section we derive exact SOET expressions for the per-site energy and double occupancy in the particular case of the uniform Hubbard system () for which the LDA decomposition of the full correlation energy in terms of per-site contributions,


is exact. Thus we extend to multiple impurities Eqs. (21) and (23) of Ref. Senjean et al. (2018).

For that purpose, let us introduce the following per-site analog of Eq. (21),


which, for a uniform density profile , gives


Note that, when combining Eqs. (22) and (23) with Eq. (21), we obtain the following expression for the bath correlation energy,


By inserting the decomposition in Eq. (24) into the exact double site-occupation expression Capelle and Campo Jr. (2013) (we denote for simplicity),


where and is the total number of electrons, it comes from Eq. (19),


By using the fact that, in the particular (uniform) case considered here, and, according to the Hellmann–Feynman theorem (see Eq. (II.1)),




for any site , it comes from the separations in Eqs. (8) and (15) that


Finally, combining Eqs. (27), (29) and (30) leads to the following exact expression for the double occupation in SOET with multiple impurities,


The expression derived in Ref. Senjean et al. (2018) in the particular case of a single impurity site is recovered from Eq. (31) when . Note also that the double occupations of the impurity sites are in principle not equal to each other in the fictitious -impurity-interacting system, simply because translation symmetry is broken, as readily seen from Eq. (12), even though the embedding potential restores uniformity in the density profile.

Turning to the per-site energy Lima et al. (2003),


where is the per-site non-interacting kinetic energy functional, we can insert Eq. (24) into Eq. (32) and use Eqs. (15) and (19), thus leading to


Moreover, applying once more the Hellmann–Feynman theorem to the variational energy expression in Eq. (II.1) gives


or, equivalently,


which, for , leads to


As a result, the second term in the right-hand side of Eq. (33) can be simplified as follows,


thus leading, according to Eq. (24), to the final exact per-site energy expression,


or, equivalently,


Note that the expression in Eq. (39), which is a generalization for multiple impurities of the energy expression derived in Ref. Senjean et al. (2018), is convenient for practical (approximate) SOET calculations where the density profile calculated self-consistently might deviate significantly from uniformity Senjean et al. (2018). Finally, as shown in Appendix B, since the exact per-site bath correlation functional fulfills the fundamental relation,


Eq. (39) can be further simplified as follows,


Let us stress that any approximate density functional of the form fulfills the exact condition in Eq. (40). This is the case for all the DFAs considered in this work. As a result, switching from Eq. (39) to Eq. (41) brings no additional errors when approximate functionals are used.

Iii Local density-functional approximations in SOET

In order to perform practical SOET calculations and compute, for example, per-site energies, we need DFAs, not only for the per-site correlation energy , like in conventional KS-DFT, but also for the per-site complementary bath correlation energy or, equivalently, for the impurity correlation energy [see Eq. (23)]. Existing approximations to the latter functionals are discussed in Secs. III.1 and  III.2, respectively. A new functional is proposed in Sec. III.3 and a simple multiple-impurity DFA is introduced in Sec. III.4. Let us stress that, in all the DFAs considered in this work, we make the approximation that the impurity correlation functional does not depend on the occupations in the bath,


or, equivalently,


where is the collection of densities on the impurity sites. The implications of such an approximation are discussed in detail in Ref. Senjean et al. (2017).

iii.1 Bethe ansatz LDA for

The BALDA approximation Lima et al. (2002, 2003); Capelle et al. (2003) to the full per-site correlation energy functional is exact for , , and all values at half-filling (). It reads as follows,


where the BALDA density-functional energy equals




The -dependent function is determined by solving

where and are zero- and first-order Bessel functions. When , the BALDA energy reduces to the (one-dimensional) non-interacting kinetic energy functional


which is exact in the thermodynamic limit ().

iii.2 Review of existing DFAs for a single impurity

iii.2.1 Impurity-BALDA

The impurity-BALDA approximation (iBALDA), which was originally formulated in Ref. Senjean et al. (2018) for a single impurity site, consists in modeling the correlation energy of the impurity-interacting system with BALDA:


In other words, the iBALDA neglects the contribution of the bath to the total per-site correlation energy [see Eq. (23) with =1],


iii.2.2 DFA based on the single-impurity Anderson model.

As shown in Ref. Senjean et al. (2018), a simple single-impurity correlation functional can be designed from the following perturbation expansion in of the symmetric single-impurity Anderson model (SIAM) correlation energy Yamada (1975),

where is the so-called impurity level width parameter of the SIAM. A density functional is obtained from Eq. (III.2.2) by introducing the following -dependent density-functional impurity level width,


A rationale for this choice is given in Ref. Senjean et al. (2018). Combining the resulting impurity correlation functional with BALDA gives the so-called SIAM-BALDA approximation Senjean et al. (2018). In summary, within SIAM-BALDA, we make the following approximations,



iii.3 Combined two-level/Bethe ansatz LDA functional

As shown in Ref. Senjean et al. (2017), in the particular case of the two-level (2L) Hubbard model (also referred to as the Hubbard dimer) with two electrons, the full density-functional correlation energy is connected to the impurity one by a simple scaling relation,


where is the occupation of the impurity site. In this case, the bath reduces to a single site with occupation . Combining Eq. (55) with BALDA gives us a new single-impurity DFA that will be referred to as 2L-BALDA in the following. In summary, within 2L-BALDA, we make the following approximations,




In our calculations, the accurate parameterization of Carrascal et al. Carrascal et al. (2015, 2016) has been used for .

iii.4 DFA for multiple impurity sites

As pointed out in Sec. II.1, in SOET, one can gradually move from KS-DFT to pure wavefunction theory by increasing the number of impurities from 0 to the number of sites. Let us, for convenience, introduce the following notation,


or, equivalently,


where is the proportion of impurity sites in the partially-interacting system. In the thermodynamic limit, becomes a continuous variable and, if the number of impurity sites is large enough, the following Taylor expansion can be used,


As readily seen from Eqs. (6), (7), (8), (20), and (25), when or, equivalently, , we have


If, for simplicity, we keep the zeroth-order term only in Eq. (60), we obtain a generalization of iBALDA for impurities, which is denoted as iBALDA() in the following. The exploration of first- and higher-order corrections in Eq. (60) is left for future work. In summary, within iBALDA(), we make the following approximations,




Iv Computational details

The various DFAs discussed in Sec. III and summarized in Table 1 have been applied to the -site uniform one-dimensional Hubbard model with an even number of electrons and . Periodic () and antiperiodic () boundary conditions have been used when [i.e. is an odd number] and [i.e. is an even number], respectively. In all the SOET calculations [note that, in the following, they will be referred to by the name of the DFA that is employed (see the first column of Table 1)], the DMRG method White (1992, 1993); Verstraete et al. (2008); Schollwöck (2011); Nakatani () has been used for solving (self-consistently or with the exact uniform density) the many-body Eq. (12). The maximum number of renormalized states (or virtual bond dimension) was set to . Standard DMRG calculations (simply referred to as DMRG in the following) have also been performed on the conventional uniform Hubbard system for comparison. For analysis purposes, exact correlation energies and their derivatives in and have been computed in a smaller 8-site ring. Technical details are given in Appendix C. The performance of the various DFAs has been evaluated by computing double occupations and per-site energies according to Eqs. (31) and (41), respectively. This required the implementation of various DFA derivatives. Details about the derivations are given in Appendices DE, and F. Note that the hopping parameter has been set to in all the calculations.

SOET method DFA used for Correlation functionals
iBALDA() Eqs. (44)–(III.1)
SIAM-BALDA Eqs. (44)–(III.1), (III.2.2) and (52)
2L-BALDA Eqs. (44)–(III.1) and (110)
Table 1: Summary of the single- and multiple-impurity Hxc DFAs used in this work for the bath. See Sec. III for further details.
Figure 1: Accurate single-impurity () correlation density-functional energies computed at the DMRG level in the half-filled case. Results obtained for are shown in black dashed lines. See text for further details.

V Results and discussion

v.1 Half-filled case

In this work, SOET is applied to a relatively small 32-site ring. Nevertheless, as illustrated in the following, the number of sites is large enough so that there are no substantial finite-size errors in the calculation of density-functional energies. This can be easily seen in the half-filled case () where the exact embedding potential in Eq. (12) equals zero in the bath and on the impurity sites (see Appendix A). The resulting impurity correlation density-functional energy [see Appendix C for further details] calculated at the DMRG level is shown in Fig. 1. The convergence towards the thermodynamic limit () is relatively fast. Note that the results obtained for and are almost undistinguishable.

Let us now focus on the calculation of double occupations which has been implemented according to Eq. (31) for various approximate functionals. Results are shown in Fig. 2.

Figure 2: Double occupations obtained from half-filled SOET calculations for various single- and multiple-impurity DFAs. Both standard reference DMRG and iBALDA(=1) results are taken from Ref. Senjean et al. (2018). See text for further details.

While, in the weakly correlated regime and up to , SIAM-BALDA gives the best results, it dramatically fails in stronger correlation regimes, as expected Senjean et al. (2018). In the particular strongly correlated half-filled case, SIAM-BALDA can be improved by interpolating between the weakly and strongly correlated regimes of the SIAM Senjean et al. (2018). Unfortunately, generalizing such an interpolation away from half-filling is not straightforward Senjean et al. (2018). On the other hand, 2L-BALDA performs relatively well for all the values of . Turning to the multiple-impurity iBALDA approximation, the accuracy increases with the number of impurities, as expected, but the convergence towards DMRG is slow. Switching from two to three impurities slightly improves on the result, which is still less accurate than the (single-impurity) 2L-BALDA one. Interestingly, the same pattern is observed in DMET when the matching criterion involves the impurity site occupation only [see the noninteracting (NI) bath formulation in Ref. Bulik et al. (2014) and the Fig. 2 therein]. While our results would be improved by designing better -dependent DFAs than iBALDA(), the performance of multiple-impurity DMET is increased when matching not only diagonal but also non-diagonal density matrix elements Knizia and Chan (2012); Bulik et al. (2014). Further connections between SOET and DMET are currently investigated and will be presented in a separate work.

Figure 3: Single-impurity density-functional correlation energy derivative with respect to calculated for various values. The derivative shown for iBALDA(=1) is, by construction, the derivative of the per-site BALDA correlation energy. Exact results were obtained for the 8-site ring. See text and Appendix C for further details.

Returning to the single-impurity DFAs (), it is quite instructive to plot the derivative in of the various functionals in order to analyze further the double occupations shown in Fig. 2. According to Eqs. (24) and (31), both full and impurity correlation energy derivatives should in principle be analyzed. However, at half-filling, and for any and values, BALDA becomes exact for when approaching the thermodynamic limit, by construction Lima et al. (2003). As a result, approximations in the impurity correlation functional will be the major source of errors which are purely functional-driven in the half-filled case, since the exact embedding potential is known (i.e. the correct uniform density profile is obtained when solving Eq. (12) in this case). Results are plotted with respect to the density in Fig. 3, thus providing a clear picture not only at half-filling but also around the latter density regime. The exact results obtained by Lieb maximization for the 8-site ring are used as reference [see the technical details in Appendix C]. Let us recall that, within iBALDA(=1), the impurity correlation energy is approximated with the full per-site correlation one , which is then modeled at the BALDA level of approximation. The substantial (negative) difference between the exact full per-site and impurity correlation energy derivatives at half-filling (see Fig. 3), which is missing in iBALDA(=1), explains why the latter approximation systematically overestimates the double occupation. Interestingly, the derivative obtained with iBALDA(=1) at , which is nothing but the derivative of the BALDA per-site correlation energy at half-filling, is essentially on top of its exact 8-site analog, thus confirming that finite-size effects are negligible. Turning to SIAM-BALDA, the derivative in of the impurity correlation energy turns out to be relatively accurate at for and , thus leading to good double occupations in this regime of correlation. The derivative deteriorates for the larger value, as expected Senjean et al. (2018). Note finally that 2L-BALDA, where the impurity correlation functional is approximated by its analog for the Hubbard dimer, is the only approximation that provides reasonable derivatives in all correlation and density regimes. The stronger the correlation is, the more accurate the method is.

Let us finally discuss the per-site energies which have been computed according to Eq. (41) for the various DFAs. Results obtained at half-filling are shown in Fig. 4.

Figure 4: Per-site energies obtained from half-filled SOET calculations for various single- and multiple-impurity DFAs. Standard DMRG results are used as reference. See text for further details.

The discussion on the performance of each DFA for the double occupation turns out to hold also for the energy. This is simply due to the fact that the per-site non-interacting kinetic energy expression in Eq. (48) is highly accurate (it becomes exact in the thermodynamic limit), like the BALDA per-site correlation energy at half-filling (and therefore its derivative with respect to ). It then becomes clear, when comparing Eqs. (31) and (41), that, at half-filling, the only source of errors in the per-site energy is, like in the double occupation, the derivative in of the impurity correlation functional.

v.2 Functional-driven errors away from half-filling

Away from half-filling, the exact embedding potential is not uniform anymore in the bath Senjean et al. (2018), thus reflecting the dependence in the bath site occupations of the impurity correlation energy or, equivalently, of the per-site bath correlation energy [see Eqs. (12), (23), and (25)]. Such a dependence is neglected in all the DFAs used in this work, which induces errors in the density when the impurity-interacting wavefunction is computed self-consistently according to Eq. (12). This generates so-called density-driven errors in the calculation of both the energy and the double occupation Kim et al. (2013). The latter are analyzed in detail in Sec. V.3.

In this section, we focus on the functional-driven errors. In other words, all density-functional contributions are calculated with the exact uniform density, like in Sec. V.1. The corresponding per-site energies are shown in Fig. 5. Only the most challenging range of fillings (i.e.  Senjean et al. (2018)) is shown.

Figure 5: Per-site energies calculated in SOET for fillings in the range with (upper curves) and (lower curves). Results obtained for various DFAs with exact (full lines) and self-consistently converged (dashed lines) densities are shown. The iBALDA(=1) and reference DMRG results are taken from Ref. Senjean et al. (2018). See text for further details.

It clearly appears that iBALDA(=1), while failing dramatically in the half-filled strongly correlated regime, performs relatively well away from half-filling in all correlation regimes. It is the best approximation in this regime of density. Increasing the number of impurities has only a substantial effect on the energy when approaching half-filling for large values. SIAM-BALDA performs reasonably well in the weakly correlated regime but not as well as the other functionals. In the same regime of correlation, 2L-BALDA stands in terms of accuracy between iBALDA and SIAM-BALDA in the lower density regime while giving the best result when approaching half-filling. The latter statement holds also in the strongly correlated regime.

In order to further analyze the performance of each functional, let us consider the derivative in of the single-impurity correlation functional shown in Fig. 3. Away from half-filling and in the strongly correlated regime, the full per-site and impurity correlation energies give the same derivative so that neglecting the last term in the right-hand side of Eq. (41), which is done in iBALDA(), is well justified. Interestingly, this feature is well reproduced by 2L-BALDA, where the impurity correlation energy obtained from the Hubbard dimer is combined with the per-site BALDA correlation energy [compare 2L-BALDA with iBALDA(=1) curves in the bottom panel of Fig. 3]. Obviously, SIAM-BALDA does not exhibit the latter feature which explains why it fails in this regime of correlation and density. Let us finally stress that, since BALDA provides an accurate description of the full per-site correlation energy in the strongly correlated regime, the second term in the right-hand side of Eq. (41) is expected to be well described in all the DFAs considered in this work (since they all use BALDA for this contribution). As shown in Fig. 6, this is actually the case [ should be compared with the iBALDA(=1) derivative].

Figure 6: Single-impurity density-functional correlation energy derivative with respect to calculated for various values. The derivative shown for iBALDA(=1) is, by construction, the derivative of the per-site BALDA correlation energy. Exact results were obtained for the 8-site ring. See text and Appendix C for further details.

Turning to the weaker correlation regime, iBALDA(=1) underestimates the derivative in of the full per-site correlation energy [compare with iBALDA(=1) in Fig. 6] away from half-filling while setting to zero the derivative in of the per-site bath correlation functional whose accurate value is actually negative [compare with the exact impurity curve in Fig. 3]. The accumulation of errors leads to the relatively accurate results shown in Fig. 5. Turning to SIAM-BALDA in the same density and correlation regime, the (negative) derivative in of the per-site bath correlation energy is significantly overestimated [compare iBALDA(=1) with SIAM-BALDA in the top panel of Fig. 3], thus giving a total per-site energy lower than iBALDA(=1), as can be seen from the upper curves in Fig. 5. Interestingly, the latter derivative in is less overestimated when using 2L-BALDA. This explains why it performs better than SIAM-BALDA but still not as well as iBALDA(=1) which benefits from error cancellations.

v.3 Self-consistent results

We discuss in this section the results obtained by solving the density-functional impurity problem self-consistently [see Eq. (12)], thus accounting for not only functional-driven but also density-driven errors Kim et al. (2013). Self-consistently converged densities obtained on the impurity site(s) for various DFAs are shown in Figs. 7 and 8.

Figure 7: Self-consistently converged densities on the impurity site(s) plotted against the filling for various DFAs and values. For symmetry reasons, the two impurities in iBALDA(=2) have the same density.
Figure 8: Self-consistently converged densities on the impurity sites plotted against the filling for iBALDA(=3) with different values. The central and neighboring impurity densities are denoted as and , respectively.

Note that, at half-filling (), the exact embedding potential has been used, thus providing, for this particular density, the exact uniform density profile. Away from half-filling, approximate density-functional embedding potentials have been used, thus giving a density profile that is not strictly uniform anymore. Interestingly, for all the DFAs except SIAM-BALDA, the deviation from uniformity is more pronounced when approaching the half-filled strongly correlated regime. In the case of iBALDA, we notice that errors in the density are attenuated when increasing the number of impurities, as expected. On the other hand, SIAM-BALDA generates huge density errors for almost all fillings when the correlation is strong.

Turning to the weaker correlation regime, we note that, in contrast to iBALDA, both self-consistent SIAM-BALDA and 2L-BALDA calculations hardly break the exact uniform density profile (see the top panel of Fig. 7). This can be rationalized by plotting the various impurity correlation potentials (see Fig. 9).

Figure 9: Single-impurity density-functional correlation potential obtained for various DFAs and values. The iBALDA(=1) potential is, by construction, the BALDA correlation potential.

Let us first stress that, for all the DFAs considered in this work, the correlation contribution to the embedding potential on site reads [see Eqs. (12), (21) (22), and (42)]


The latter potential will therefore be equal to zero on the impurity sites if iBALDA is used. It will be uniform and equal to (at least in the first iteration of the self-consistent procedure as we start with a uniform density profile) in the bath. As shown in the top panel of Fig. 9 [see the iBALDA(=1) curve], in the weakly correlated regime, the latter potential, which is described with BALDA, is strongly attractive for densities lower than 0.6. This will induce a depletion of the density on the impurity sites, as clearly shown in the top panels of Figs. 7 and 8. The opposite situation is observed for densities in the range , as expected from the strongly repulsive character of the potential. Note that this charge transfer process, which is completely unphysical Akande and Sanvito (2010), is due to the incorrect linear behavior in of the BALDA correlation potential away from the strongly correlated regime Senjean et al. (2018). As readily seen from the first-order expansion in given in Eq. (32) of Ref. Senjean et al. (2018), the latter potential is expected to change sign at , which is in complete agreement with the top panel of Fig. 9. Returning to Eq. (64), for the single-impurity DFAs SIAM-BALDA and 2L-BALDA, the correlation contribution to the embedding potential reads on the impurity site and still in the bath. As shown in the top panel of Fig. 9, SIAM-BALDA and 2L-BALDA impurity correlation potentials dot not deviate significantly from zero and, unlike iBALDA, they do not exhibit unphysical features in the weakly correlated regime, which explains why both DFAs give relatively good densities. Turning to iBALDA(=3) in the strongly correlated regime (middle and bottom panels of Fig. 8) and in the vicinity of half-filling, the impurity site 1 better reproduces the physical occupation than its nearest impurity neighbors (sites 0 and 2). It can be explained as follows. Site 1 is the central site of the (=3)-impurity-interacting fragment, while sites 0 and 2 are directly connected to the (non-interacting and therefore unphysical) bath. As a consequence, site 1 “feels” the bath less than sites 0 and 2. It is, like in the physical system, surrounded by interacting sites. Interestingly, a similar observation has been made by Welborn et al. Welborn et al. (2016) within the Bootstrap-Embedding method, which is an improvement of DMET regarding the interaction between the fragment edge and the bath.

Let us now refocus on the poor performance of SIAM-BALDA in the strongly correlated regime. As shown in the middle and bottom panels of Fig. 9, the correlation part of the embedding potential on the impurity site [which corresponds to the difference between iBALDA(=1) and SIAM-BALDA curves] is, in this case, repulsive for densities lower than about 0.25 and strongly attractive in the range . The latter observation explains the large increase of density on the impurity site in the self-consistent procedure, as depicted in the bottom panel of Fig. 7.

Finally, as shown in Fig. 5, using self-consistently converged densities rather than exact (uniform) ones can induce substantial density-driven errors on the per-site energy, especially in the strongly correlated regime. Interestingly, both functional- and density-driven errors somehow compensate for 2L-BALDA around half-filling (see the lower curves in Fig. 5). This is not the case anymore exactly at half-filling since we used the exact embedding potential, thus removing density-driven errors completely. Regarding SIAM-BALDA, the large error in the converged density obtained for at (see the middle panel of Fig. 7) is reflected on the per-site energy. In this case, errors just accumulate. For larger fillings, SIAM-BALDA gives a better density on the impurity site but the functional-driven errors remain substantial.

v.4 Derivative discontinuity at half-filling

As illustrated in Ref. Senjean et al. (2018) in the special case of the atomic limit (i.e. or, equivalently, ), the impurity correlation potential on the impurity site should undergo a discontinuity at half-filling in the thermodynamic limit (even for finite values of ). This feature has fundamental and practical implications, in particular for the calculation of the fundamental gap. The latter quantity plays a crucial role in the description of the metal-insulator transition which, in the one-dimensional Hubbard model, appears as soon as the on-site repulsion is switched on (Lieb and Wu (1968).

As shown in Fig. 9, it is present in the iBALDA and SIAM-BALDA functionals, as a consequence of the hole-particle symmetry condition used in their construction. For the former functional, the feature is simply inherited from BALDA. On the other hand, even though the 2L-BALDA impurity correlation potential, which is extracted from the Hubbard dimer, also fulfills the particle-hole symmetry condition, it smoothly tends to 0 when approaching . As a consequence, the potential exhibits no derivative discontinuity (DD) at half-filling when is finite. It only does when  Senjean et al. (2017). As discussed by Dimitrov et al. Dimitrov et al. (2016), this is due to the fact that the dimer functional reproduces an intra-system steepening and not an inter-system DD. In other words, the change in density in the functional does not correspond to a change in the total number of electrons. The latter is indeed fixed to 2 in the dimer. Only the number of electrons on the impurity site varies. The problem becomes equivalent to describing an inter-system DD only when the impurity can be treated as an isolated system, which is indeed the case in the atomic (or, equivalently, strongly correlated) limit.

Note finally that, from a practical point of view, exhibiting a DD is not necessarily an advantage as convergence problems may occur around half-filling Lima et al. (2003). In this work, this problem has been bypassed by using the exact embedding potential at half-filling. Also, given that only 32 sites are considered, the closest uniform occupation to half-filling is obtained for 30 electrons, i.e. , which is far enough from the strictly half-filled situation. Convergence issues are expected to arise when approaching the thermodynamic limit since the density can then be much closer to 1. The practical solutions to this problem, which have been proposed for conventional KS calculations and use either a finite temperature Xianlong et al. (2012) or ad-hoc parameters Kurth et al. (2010); Karlsson et al. (2011); Ying et al. (2014), could in principle be implemented in SOET. This is left for future work. Note finally that, despite the absence of DD in the 2L-BALDA impurity correlation potential, the BALDA correlation potential is still employed for the bath within 2L-BALDA, which means that convergence issues will appear as soon as occupations in the bath are close to 1.

Vi Conclusions and perspectives

Several extensions of a recent work Senjean et al. (2018) on site-occupation embedding theory (SOET) for model Hamiltonians have been explored. Exact expressions for per-site energies and double occupations have been derived for an arbitrary number of impurity sites. A simple -impurity embedding density-functional approximation (DFA) based on the Bethe ansatz local density approximation (BALDA) and referred to as iBALDA() has been proposed and tested on the one-dimensional Hubbard Hamiltonian. A new single-impurity DFA [referred to as 2L-BALDA] which combines BALDA with the impurity correlation functional of the two-level (2L) Hubbard system Senjean et al. (2017) has also been proposed. Finally, the performance of an existing DFA based on the single impurity Anderson model (SIAM) and BALDA Senjean et al. (2018), hence its name SIAM-BALDA, has been analyzed in further details in all correlation regimes. Both functional- and density-driven errors have been scrutinized. Among all the single-impurity DFAs, 2L-BALDA is clearly the one that performs the best in all density and correlation regimes. Unfortunately, the convergence of the (too) simple iBALDA approximation in the number of impurities towards the accurate DMRG results was shown to be slow. Better DFAs for multiple impurities are clearly needed. This is left for future work.

SOET can be extended further in many directions. First of all, substituting a Green function calculation for a many-body wavefunction one like DMRG is expected to reduce the computational cost of the method. From a formal point of view, this would also enable us to connect SOET with the dynamical mean-field theory (DMFT) [see, for example, Ref. Ayral et al. (2017)]. Note that the current formulation of SOET is canonical. It would be interesting to remap (density wise) the original fully-interacting Hubbard problem onto an open impurity system, in the spirit of the density matrix embedding theory (DMET). This may be achieved, in principle exactly, by combining SOET with partition DFT Elliott et al. (2010), or, in a more approximate way, by projecting the whole impurity-interacting problem in SOET onto a smaller embedded subspace. A Schmidt decomposition could be employed in the latter case Knizia and Chan (2012). Finally, an important step forward would be the exploration of SOET in higher dimensions. Work is currently in progress in all these directions.

Let us finally mention that the basic idea underlying SOET, which consists in extracting site (or orbital) occupations from a partially-interacting system, can be extended to an ab initio quantum chemical Hamiltonian by using, for example, its simplified seniority-zero expression Richardson (1963); Richardson and Sherman (1964); Limacher (2016). This will be presented in a separate work.


E.F. would like to dedicate this work to the memory of János G. Ángyán. He would also like to thank Andreas Savin for a fruitful discussion on the train from Middelfart to Copenhagen. B.S. thanks D. Carrascal for taking the time to check the parameterization of his Hubbard dimer functional, and M. Saubanère, L. Mazouin, and K. Deur for fruitful discussions. This work was funded by the Ecole Doctorale des Sciences Chimiques 222 (Strasbourg), the ANR (MCFUNEX project, Grant No. ANR-14-CE06- 0014-01), the “Japon-Unistra” network as well as the Building of Consortia for the Development of Human Resources in Science and Technology, MEXT, Japan for travel funding.

Appendix A Exact embedding potential at half-filling for multiple impurities

Let us consider any density summing up to a number of electrons. Under hole-particle symmetry, this density becomes and the number of electrons equals where is the number of sites. We will prove that these two densities give the same correlation energy for the -impurity interacting system. Since, for any local potential , the variational principle in Eq. (5) reads as follows for an impurity-interacting system,


which gives, for any density ,


thus leading to the Legendre–Fenchel transform expression,


By applying a hole-particle symmetry transformation to Eq. (67) [we will now indicate the number of particles in the impurity-interacting energies for clarity], we obtain

where is the ()-particle ground-state of the following -impurity interacting Hamiltonian:


Applying the hole-particle transformation to the creation and annihilation operators,


to the -impurity-interacting Hamiltonian in Eq. (69) leads to


or, equivalently,

Then, by substituting and shifting the potential as follows,


we finally obtain