Multiple impurities and combined local density approximations in Site-Occupation Embedding Theory
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.
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
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,
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),
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)
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
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.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,
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,
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,
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
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)),
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),
Moreover, applying once more the Hellmann–Feynman theorem to the variational energy expression in Eq. (II.1) gives
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,
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,
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
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
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,
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,
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,
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 D, E, and F. Note that the hopping parameter has been set to in all the calculations.
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
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.
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
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
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
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
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.
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
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.
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].
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.
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).
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
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.
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
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
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
Then, by substituting and shifting the potential as follows,
we finally obtain