Exchange-correlation potentials for multi-orbital quantum dots subject to generic density-density interactions and Hund's rule coupling
EExchange-correlation potentials for multi-orbital quantum dots subject to genericdensity-density interactions and Hund’s rule coupling
Nahual Sobrino,
1, 2
Stefan Kurth,
2, 3, 1 and David Jacob
2, 3, ∗ Donostia International Physics Center (DIPC),Paseo Manuel de Lardizabal 4, E-20018 San Sebasti´an, Spain Nano-Bio Spectroscopy Group and European Theoretical Spectroscopy Facility (ETSF),Dpto. de F´ısica de Materiales, Universidad del Pa´ıs Vasco UPV/EHU, Av. Tolosa 72, E-20018 San Sebasti´an, Spain IKERBASQUE, Basque Foundation for Science, Maria Diaz de Haro 3, E-48013 Bilbao, Spain (Dated: August 10, 2020)By reverse-engineering from exact solutions we obtain Hartree-exchange-correlation (Hxc) po-tentials for a double quantum dot subject to generic density-density interactions and Hund’s rulecoupling. We find ubiquitous step structures of the Hxc potentials that can be understood andderived from an analysis of stability diagrams. We further show that a generic Hxc potential canbe decomposed into four basic potentials which allows for a straight-forward parametrization andpaves the road for the construction of Hxc potentials for interacting multi-orbital systems. Finallywe employ our parametrization of the Hxc potential in density functional theory calculations ofmulti-orbital quantum dots and find excellent agreement with exact many-body calculations.
I. INTRODUCTION
Density functional theory (DFT) is one of the mostsuccessful and popular approaches for computing theelectronic structure of molecules and solids.
Its suc-cess is largely owed to its relative simplicity as well as itslow computational cost as compared to other quantummany-body approaches. While DFT is an in principleexact theory for computing the ground state energy anddensity of a many-electron system, in practice approx-imations have to be made to the exchange-correlation(xc) part of the total energy functional. The most pop-ular approximations are the local density (LDA) andgeneralized gradient approximations (GGA) in con-densed matter physics, and the so-called hybrid func-tionals in chemistry. While these approximations usu-ally work quite well for systems with weak to moderateelectronic interactions, they completely fail for so-calledstrongly correlated systems where the interactions be-tween electrons dominate over the kinetic energy. SinceDFT is a formally exact theory which is valid also in thestrongly correlated regime, this failure has to be assignedto shortcomings of the approximations used. Apparently,the crucial ingredient missing in standard functionals isthe so-called derivative discontinuity, i.e., the discon-tinuous jump of the exact xc potential of an open sys-tem as the particle number crosses an integer. Instrongly correlated systems the derivative discontinuitycontributes a substantial part, e.g., to the fundamentalgap or plays a crucial role in the binding and dissocia-tion of molecules. The essential physics of strongly correlated systemscan already be captured by simple but highly nontriviallattice models such as the Anderson or the Hubbardmodel which can be solved by advanced many-bodymethods as for example Dynamical Mean-Field Theory(DMFT). One way around the problems of the stan-dard approximations of DFT is then to combine DFT with advanced many-body calculations of lattice models.In these approaches the model Hamiltonian describes thestrongly interacting part of the system, while the weaklyto moderately interacting part is still described at thelevel of DFT. An important example is the combinationof DFT with DMFT (DFT+DMFT). Originally estab-lished for the description of bulk materials, more recentlythe DFT+DMFT approach has been applied to the de-scription of nanoscale systems and molecules . How-ever, this approach is hampered by the so-called double-counting problem , limiting its predictivity. More re-cently, however, new efforts in combining DFT with lat-tice models avoiding the double-counting problem orto solve the double-counting problem in DFT+DMFT have been undertaken.Lattice models can also be solved by a lattice ver-sion of DFT, an idea which has been pioneered bySch¨onhammer. Later this approach has been fur-ther extended to study lattice problems not only inequilibrium but also in out-of-equilibrium situationssuch as external time-dependent driving fields or(steady-state) transport . A common theme in manyof these studies is again the crucial role of the deriva-tive discontinuity (or, alternatively, step features in thexc potential) in the correct description of strongly corre-lated systems. For instance, for the Hubbard model thederivative discontinuity is the only contribution to theband gap, i.e. the mechanism responsible for openingthe Mott-Hubbard gap within DFT. Here, in the same framework of lattice DFT, we aim fora better understanding of the structure of the (equilib-rium) xc potentials of small multi-orbital models. This ismotivated by earlier studies on the single-impurity An-derson model , i.e. a single interacting impurity incontact with non-interacting leads. There it was shownthat the Kondo effect, one of the hallmarks of strong cor-relation in the context of electron transport, can be de-scribed at the DFT level if the corresponding xc potentialexhibits a step at integer particle number on the impu- a r X i v : . [ c ond - m a t . m e s - h a ll ] A ug rity. This step is already present in the xc potential of anisolated impurity in contact with a particle and heat bathdescribed via the grand canonical ensemble . The exactxc potential of the impurity connected to leads has thesame qualitative features as the uncontacted one with thequantitative difference that the broadening of the step isdetermined (at low temperature) by the coupling to theleads while in the uncontacted case it is determined bythe temperature of the bath. Therefore, we expect thatknowledge of the xc potential of multi-orbital models inequilibrium, apart from being interesting in itself, have adirect relevance in the context of transport.To this end we extensively study a double quantum dot(DQD) subject to generic interactions as the simplestpossible model system for a strongly correlated multi-orbital system. Using reverse-engineering, we constructthe exact xc potentials whose essential features are stepstructures which depend on the particular choice of theinteraction. We illustrate how these step structures canbe inferred from an analysis of the stability diagram, i.e.,regions of different ground states in the parameter spacegiven by the single-particle level structure. Analysis ofthe reverse-engineered xc potentials reveals that they canbe constructed from four basic building blocks which canbe rationalized by a corresponding decomposition of theelectron-electron interaction. It is exactly this decom-position of the interaction and the corresponding xc po-tentials which then allows us to build xc functionals formulti-orbital quantum dots with more than two orbitals. II. MODEL
We consider a multi-orbital quantum dot (QD) with M orbitals subject to direct Coulomb repulsion and Hund’srule coupling. The corresponding Hamiltonian reads: H = (cid:88) α v α ˆ n α + (cid:88) α U α ˆ n α ↑ ˆ n α ↓ + (cid:88) α<β U αβ ˆ n α ˆ n β − (cid:88) α<β,σ J αβ (cid:104) ˆ n ασ ˆ n βσ + (cid:16) c † ασ c α ¯ σ c † β ¯ σ c βσ (cid:17)(cid:105) (1)where c ασ ( c † ασ ) are the annihilation (creation) opera-tors for orbital α and spin σ , ˆ n ασ is the correspond-ing number operator and ˆ n α = ˆ n α ↑ + ˆ n α ↓ . U α isthe direct intra-orbital Coulomb repulsion for orbital α , U α ≡ (cid:104) α, α | ˆ V c | α, α (cid:105) , and U αβ is the direct inter-orbitalCoulomb repulsion between electrons in two different or-bitals, U αβ = (cid:104) α, β | ˆ V c | α, β (cid:105) for α (cid:54) = β . J αβ is the Hund’srule coupling, i.e. the exchange integral of the Coulombinteraction, J αβ = (cid:104) α, β | ˆ V c | β, α (cid:105) for α (cid:54) = β . Note thathere we have split already the Hund’s rule term into thedensity-density contribution (first term in the last line)and the spin-flip contrbitution (last term). v α are theorbital energies (single-particle energies) of the orbitals α which can be tuned by an external “gate” potential.From here on we will thus refer to v α as the gate poten-tial or simply gate for orbital α . Here we work at (typically small) finite temperature T and consider the grand canonical ensemble (GCE).The corresponding density matrix (statistical operator)is thus given byˆΓ = e − β H Z = 1 Z (cid:88) m e − βE m | m (cid:105) (cid:104) m | (2)where β = 1 /T and Z is the GCE partition function withthe chemical potential µ set to zero for convenience. The | m (cid:105) are the many-body eigenstates of the QD and E m thecorresponding eigenenergies, i.e. H | m (cid:105) = E m | m (cid:105) . Notethat in the absence of Hund’s rule coupling ( J αβ = 0)the many-body eigenstates | m (cid:105) are simply Slater deter-minants built from the single particle orbitals | φ ασ (cid:105) = c † ασ | (cid:105) where | (cid:105) is the vacuum state.The Hamiltonian (1) is very common in the fields ofstrongly correlated electrons and mesoscopic physics, asit provides a natural description of 3d- or 4f-shells of tran-sition metal or lanthanide impurities in metallic hostsand of multi-orbital quantum dots. In these systemsdensity-density interactions and Hund’s rule coupling areby far the most important interactions. In particular therole of the latter has become a focus of intense research inthe field of strongly correlated electrons in recent years. Alternatively, the Hamiltonian (1) may be viewed as alattice Hamiltonian in the limit of vanishing hopping be-tween sites which is similar to the point of view taken instrictly correlated DFT.
III. REVERSE-ENGINEERING OF HXCPOTENTIALS AT FINITE TEMPERATURE
In our truncated Hilbert space the density is uniquelydefined by all occupancies n ≡ ( n , . . . , n M ) of the QDorbitals. The Mermin theorem (the finite-temperatureversion of the Hohenberg-Kohn theorem ) then estab-lishes a one-to-one correspondence between the den-sity (occupancies) n and the external potential (gate) v ≡ ( v , . . . , v M ): n − ←→ v or, in other words, the ex-ternal potential is a functional of the occupancies, i.e., v α = v α [ n ].In order to proceed we introduce the Kohn-Sham (KS)system i.e. an effective non-interacting system that ex-actly reproduces the density n of the many-body Hamil-tonian H . Here the KS Hamiltonian is already diagonalin the original single-particle basis, i.e., H s = (cid:88) α v sα ˆ n α . (3)and the KS orbitals are identical to the original basisorbitals | φ sασ (cid:105) ≡ | φ ασ (cid:105) with their eigenenergies givenby the KS (gate) potentials v sα . The Hartree-exchange-correlation (Hxc) potentials v Hxc α are defined as the differ-ence between the KS gate and the actual gate potential: v Hxc α [ n ] = v sα [ n ] − v α [ n ] . (4)The Hxc potential depends on the electron density whichis completley determined by the occupancies of the QDorbitals n .In order to determine the Hxc potential v Hxc as a func-tional of the density n , the many-body problem given by H is solved for a given set of gates v . The resulting setof eigenstates and corresponding energies determines thedensity in the GCE according to: n α = Tr[ˆΓ ˆ n α ] = 1 Z (cid:88) m (cid:104) m | ˆ n α | m (cid:105) e − βE m . (5)The density in turn uniquely determines the KS potentialand thus the Hxc potential. In our case of an isolatedQD at finite temperature the occupancy n α is simplydetermined by the gate v sα of a non-interacting QD, andis thus simply given by the Fermi-Dirac distribution, i.e. n α = 2 f ( v sα ). Hence the KS gate for orbital α is given by v sα = β ln (cid:16) n α − (cid:17) and the corresponding Hxc potentialcan be obtained using (4) as: v Hxc α = 1 β ln (cid:18) n α − (cid:19) − v α . (6)Hence we have found the mapping n −→ v Hxc . By ex-ploring the parameter space v = ( v , . . . , v M ) we canestablish this mapping for the entire space of densities n = ( n , . . . , n M ) (for n α ∈ [0 , A. Hxc potentials and link to stability diagramsfor the double quantum dot
We now focus on the two-orbital case, i.e. a doublequantum dot (DQD) with generic density-density inter-actions ( U , U , U ). For the moment we neglect Hund’srule coupling ( J H = 0), but we will discuss the effect of fi-nite J H later in Sec. IV C. The Hxc potentials can be con-structed by reverse engineering as explained above. Herewe are interested in the qualitative structure of the Hxcpotentials, in particular in the positions (and heights)of step structures which appear in the low-temperaturelimit. In fact, these steps are not only the crucial but alsothe only features of the Hxc potential in the limit of lowtemperatures. In this section we will show how these stepstructures can be deduced completely from the stabilitydiagrams.A stability diagram highlights the occupations (densi-ties) of the ground states in the different regions of theplane of external gates v and v . The position and shapeof these regions in the v - v plane depend on the valuesof the interaction parameters but within each region thepair of densities ( n , n ) remains constant at (close to)zero temperature and the possible values of the local den-sities n i are restricted to 0, 1, and 2. For general tem-peratures the domain of physically realizable densities isrestricted to the square 0 ≤ n i ≤
2. In this domain ofdensities, each of the nine pairs of densities ( n , n ) with -4-2 0 -4 -2 0 v v (a) (0,0)(1,0)(2,0)(2,1) (0,1)(0,2)(1,2)(2,2)
0 1 2 n (b) v Hxc,1
01 0 1 2 n n v Hxc,2 -6-3 0 -6 -3 0 (d)v (0,0)(1,0)(2,0)(2,1) (0,1)(0,2)(1,1)(1,2)(2,2)
0 1 2 (e) v Hxc,1 U U +U U +2U
01 0 1 2(f)n v Hxc,2 U +U U U +2U FIG. 1. Panels (a)-(c) (constant interaction model, CIM):stability diagram (a) and Hxc potentials for orbitals 1 and2 (panels (b) and (c), respectively) of the double quantumdot for U = U = U . Panels (d)-(f) (Regime I): stabilitydiagram (d) and Hxc potentials for orbitals 1 and 2 (pan-els (e) and (f), respectively) of the double quantum dot for U = 2 . U , U = 3 U . All energies in units of smallestinteraction ( U ). n i ∈ { , , } corresponds to a single point which we calla vertex. Therefore a region of constant density in thestability diagram (i.e., in the v - v plane) directly cor-responds to a vertex (a single point) in the domain ofrealizable densities. A similar duality between regionsin the plane of “potentials” and vertices in the plane of“densities” has also been observed in the framework ofsteady-state DFT (i-DFT). It turns out that the structure of the Hxc potentials(in the limit of low temperatures) for a given set of inter-action parameters can be extracted just by looking at thestability diagram: for a given pair of ground state den-sities (or vertex in the density domain) one just needsto find all adjacent regions corresponding to a differentvertex. Then the Hxc potentials, which are functions ofthe density, will only contain steps which connect a givenvertex (in the density domain) with those vertices cor-responding to directly adjacent regions. The heights ofthese steps can also be extracted from the stability dia-gram. Below we will illustrate how this works presentingsome representative examples and we will also explainthe physical reasons behind our observations.As a first example we choose a simple one where allinteraction parameters are equal, U = U = U . In thiscase the total interaction can be written as U ˆ N ( ˆ N − N = ˆ n + ˆ n is the operator for the total numberof electrons on the dot. This model is known as the con-stant interaction model (CIM). It can be shown thatat zero temperature the Hxc potential v Hxc α of the CIMis independent of α and is a piecewise constant functionof the total electron number N with discontinuous stepsof height U whenever N crosses an integer. We mentionthat the CIM Hxc potential is strictly discontinuous onlyat zero temperature (this is a manifestation of the famousderivative discontinuity of DFT ). At finite but smalltemperature, the step structure persists but the Hxc po-tentials are now continuous functions of the densities .We now show how the known CIM Hxc potentials (atlow temperature) can be inferred directly from the sta-bility diagram. This diagram is shown in panel (a) ofFig. 1 for the CIM with U = U = U = 1. Herethe regions corresponding to the different possible groundstate densities (given in parenthesis) are marked by dif-ferent colors. The reverse-engineered Hxc potentials fororbitals 1 and 2 are shown in panels (b) and (c) of Fig. 1,respectively. In the stability diagram, the domain cor-responding to the occupation (0 ,
0) is directly adjacentonly to the domains with occupations (1 ,
0) and (0 , ,
0) vertex with one of those vertices inthe n - n plane we see that the resulting lines run alongthe border of the allowed density domain. The completeset of lines along the borders of the density domain fol-low from the sequence of vertices (0 , → (1 , → (2 , , → (0 , → (0 , , → (2 , → (2 , , → (1 , → (2 , v − v plane are (i) (1 , → (0 , , → (0 , , → (1 , N = n + n in the Hxcpotentials, see panels (b) and (c) of Fig. 1. Moreover,the height of these steps can also be deduced from thestability diagram: consider a v > v . Then we haveessentially a single-site model (SSM) in contact with aparticle and heat bath because the second dot doesn’tcontribute. However, we know that the the Hxc potentialof a SSM in the low temperature limit is a step functionwith step of height U at half filling . In a self-consistentDFT calculation this step in the Hxc potential leads to apinning of the KS level to the Fermi energy over a rangeof gates of width U . Therefore, the width (in v ) of the(1 ,
0) region is just U = U . Similarly, the width (in v )of the (0 ,
1) region is U = U . The line in the stabilitydiagram where the (1 ,
0) and (0 ,
1) regions touch is theline v = v for which the states with the correspondingoccupations are degenerate. The KS system is a systemof effectively non-interacting electrons which reproducesthe interacting density. However, for a non-interactingdouble dot with potentials v s and v s the only possibilityfor the half-filled dots to be degenerate is for the point v s = v s = ln 3 /β which in the limit of zero temperature approaches v s = v s = 0. Therefore, in order to reproducethe degeneracy as observed in the stability diagram, theKS potential for both orbitals has to be pinned over aninterval of range U . This can only be achieved if the Hxcpotentials for both orbitals have steps at N = 1 of height U as observed in the reverse-engineered Hxc potentials.On the line v = v also the states with occupations(2 ,
0) and (0 ,
2) as well as the ones with occupations (2 , ,
2) are degenerate. It is easy to show that alongthis line the states with occupations (2 ,
0) and (0 ,
2) arelowest in energy for the region − U > v > − U . Fornon-interacting systems, again there is only one point( v s = v s = 0) for which the states with (2 ,
0) and (0 , − U > v > − U ,the KS potentials are pinned to zero. This can only beachieved if both Hxc potentials exhibit a step of height U at N = 2, as observed. Finally, there is yet another stepof height U in both Hxc potentials for N = 3 which fol-lows in a similar way from the analysis of the contact linebetween the (2 ,
1) and (1 ,
2) regions. In this way we havetherefore been able to reconstruct the Hxc potentials ofthe CIM only by analyzing the stability diagram.In the second example we make all the interaction pa-rameters different from each other, i.e., the levels arenow not equivalent any more. Moreover, we choose theinterdot interaction U to be smaller than both U and U . This parameter regime ( U < U , U ) we denote asRegime I, see discussion in Section IV A. In the stabilitydiagram for this regime (panel (d) of Fig. 1) we now havea new region with densities (1 ,
1) showing up. In orderto deduce the low-temperature Hxc potentials (reverse-engineered results shown in panels (e) and (f) of Fig. 1),we begin by looking at the regions with occupations (1 , , v = v and for − U < v < v = v inthe same interval, we need the KS potentials on both or-bitals to be pinned to the Fermi energy. Therefore bothHxc potentials need to exhibit a step of height U alongthe line connecting the vertices (1 ,
0) and (0 , v Hxc1 ( n ,
0) = U for 1 < n < v Hxc2 (0 , n ) = U for 1 < n <
2. The regions (1 ,
0) and (1 ,
1) are adja-cent along the line v = − U for − U < v < − U and thus the KS potential of the first orbital needs to bepinned to the Fermi energy for this range of v leadingto a step of height U − U along the line connecting the(1 ,
0) and (1 ,
1) vertices for v Hxc1 . Similarly, v Hxc2 needsto exhibit a step of height U − U along the line con-necting the (0 ,
1) and (1 ,
1) vertices. Next, the regions(2 ,
0) and (1 ,
1) are adjacent for − U − U α < v α < − U α ( α = 1 ,
2) and therefore both Hxc potentials have a stepof height U along the lines connecting the (2 ,
0) and(1 ,
1) vertices. Similarly, there also has to be a step of -8-4 0 -8 -4 0 v v (a) (0,0)(1,0)(2,0)(2,1) (1,1) (0,1)(0,2)(1,2)(2,2)
0 1 2 n (b) v Hxc,1 U U +U U +2U
01 0 1 2 n n v Hxc,2 -U U U U +U +U -U U +2U -8-4 0 -8 -4 0 (d)v (0,0)(1,0)(2,0)(2,1) (0,1)(0,2)(1,2)(2,2)
0 1 2 (e) v Hxc,1 +U )/22U U +(U -U )/2U +2U
01 0 1 2(f)n v Hxc,2 +U -U U U +2U FIG. 2. Panels (a)-(c) (Regime II): stability diagram (a) andHxc potentials for orbitals 1 and 2 (panels (b) and (c), re-spectively) of the double quantum dot for U = 4 U and U = 2 U . Panels (d)-(f) (Regime III): stability diagram (d)and Hxc potentials for orbitals 1 and 2 (panels (e) and (f),respectively) of the double quantum dot for U = 2 U , and U = 2 . U . All energies in units of the smallest interaction( U ). height U in both Hxc potentials along the line con-necting the (0 ,
2) and (1 ,
1) vertices. The regions (1 , ,
2) are adjacent for − U < v < − U − U leading to a step of height U − U in v Hxc1 along theline (1 , → (1 , U − U in v Hxc2 along the (1 , → (2 ,
1) line. Finally,the regions (2 ,
1) and (1 ,
2) are adjacent along a line oflength U leading to a step of this height in both Hxcpotentials along the (2 , → (1 ,
2) line. In this way wenow have completely determined the (low temperature)Hxc potentials of both orbitals just by analyzing the sta-bility diagram. Their overall structure is such that theyexhibit steps for integer total occupation N = n + n forboth Hxc potentials plus an additional step at n α = 1 for v Hxc α . Note also that for the special case U = 0 only thesteps at n α = 1 for v Hxc α survive while those at integer N disappear. This is not surprising since in this case ourmodel just describes two completely independent singleimpurities and, naturally, the corresponding Hxc poten-tial for orbital α is completeley independent of the otherorbital and given by the Hxc potential of a SSM withinteraction strength U α . This has also been discussed as“intra-system steps” in Ref. 52. We have identified two further parameter regimes forthe interaction parameters (see Section IV A) where qual-itative changes both in the stability diagram as wellas in the Hxc potentials occur. In both regimes theinter-orbital interaction U is smaller than at leastone of the intra-orbital ones. Without loss of general-ity we may assume U ≤ U . In Regime II we have U < U < ( U + U ) / U ≤ ( U + U ) / ≤ U . In panels (a)-(c) of Fig. 2we show the stability diagram and Hxc potentials for in-teraction parameters chosen in Regime II. Compared toRegime I (panels (d)-(f) of Fig. 1), in the stability dia-gram we now find that there exists a range of potentialsfor which regions (2 ,
0) and (0 ,
1) are directly adjacentand, similarly, for the regions (2 ,
1) and (0 , n = 1 (present in Regime I)now disappears while in v Hxc , the step at n = 1 sur-vives (this step is related to the vertical lines delimitingthe (1 ,
1) region in the stability diagram). We have an-notated the plateau values in both Hxc potentials whichcan be found by analyzing the stability diagram usingsimilar arguments to the ones used above for Regime I.Finally, in panels (d)-(f) of Fig. 2 we show the stabilitydiagram and Hxc potentials for interaction parameterschosen in Regime III. Compared to Regime II, the mainqualitative difference is the disappearance of the step at n = 1 in the Hxc potential of orbital 2. Again, all thestep structures in the Hxc potentials can fully be deducedby analyzing the stability diagram.Before we close this section we would like to mentionthat for multilevel dots beyond the double dot studiedhere, in principle the analysis of the (multidimensional)stability diagram(s) also allows for a complete deductionof the low-temperature Hxc potentials of the differentorbitals. However, it is clear that this procedure rapidlybecomes quite cumbersome as the number of levels in-creases. IV. MODELLING OF THE HXC POTENTIALSA. Decomposition of the interaction into basicbuilding blocks
In the following we show that the Hxc potentials forgeneric density-density interactions can be built froma few basic potentials. We start with the most com-mon (or natural) situation where the inter-orbital inter-action U is smaller than both of the intra-orbital ones, U ≤ U , U . A specific case with U < U < U wasstudied in the previous section [see Fig. 1(d-f)]. The cor-responding Hxc potential shows steps at integer valuesof N = n + n , connected to a CIM potential, as wellas steps at n = 1 for orbital 1 or at n = 1 for orbital2 connected to a SSM potential of the corresponding or-bital.This suggests that in the regime U ≤ U , U the Hxcpotential for each orbital may be built from a superpo-sition of a CIM potential plus a SSM potential. We canrationalize this idea by a decomposition of the Coulombinteraction term as follows. Rewriting the inter-orbitalrepulsion as U ˆ n ˆ n = U N ( ˆ N − − U (cid:88) α ˆ n α ↑ ˆ n α ↓ (7)we can split the interaction into a CIM part and two SSMinteractions (one for each orbital): V int = 12 U ˆ N ( ˆ N −
1) + (cid:88) α δU α ˆ n α ↑ ˆ n α ↓ (8)where δU α ≡ U α − U is the “excess interaction” for eachorbital. This suggests to write the Hxc potential for level α for Regime I ( U ≤ U , U ) as the sum of the CIMHxc potential for interaction U and the SSM potentialfor δU α : v Hxc α [ n ] = v HxcCIM ( U )[ N ] + v HxcSSM ( δU α )[ n α ] (9)where N = n + n .Now if U is larger than at least one of the intra-orbitalinteractions U α this decomposition of the Coulomb in-teraction obviously leads to negative interactions δU α inthe SSM parts. Since the step in the Hxc potential ofthe SSM at n α = 1 would actually vanish for negativeinteractions, in this regime the step structure can ob-viously not be rationalized by the above decompositionof the interaction. Indeed the structure of the reverse-engineered Hxc potentials (Fig. 2) appears to be quitedifferent from that for the regime U ≤ U α . Essentially,two new features are found in this regime: (i) an increaseof the step height at N = 2 with respect to the CIM po-tential, and (ii) peculiar new steps at integer values of n / n . The steps at integer n / n are generatedby a peculiar interaction of the form V skew = U n ( ˆ N −
1) (10)which we will refer to as
Skew interaction from now on.This interaction is realized by setting U = 0 and U = U / U/
2, The corresponding stability diagram andthe Hxc potential for orbital 2 is shown in Fig. 3(b). Notethat the Hxc potential of orbital 1 has the same structurebut the step heights are lower by a factor of 1/2.Common to all cases is that there is always a contri-bution of the CIM potential, as long as all interactions( U , U , U ) remain finite. This contribution is given bythe smallest interaction. In the case that U is largerthan at least one of the intra-orbital interactions U α , wemay assume without loss of generality that U is thesmallest interaction. Subtracting the CIM interaction -6-4-2 0 2 4 -4 -2 0 2 4(a) v v (0,0)(0,1)(0,2)(2,0)(2,1)(2,2) 0 1 2 0 1 2(b) n n -6-4-2 0 2 4 -6 -4 -2 0 2 4(c) v v (0,0)(2,0) (0,2)(2,2) 0 1 2 0 1 2(d) n n FIG. 3. (a,b) Stability diagram (a) and Hxc potential oforbital 2 (b) for the Skew interaction V skew = U ˆ n ( ˆ N − U/ U ). (c,d) Stability diagram (c) and Hxc potential of bothorbitals (d) for the inter-orbital interaction V inter = U ˆ n ˆ n .All energies in units of U in both cases. ∼ U from the total interaction thus yields: V int − U N ( ˆ N −
1) == ( U − U ) n n + ( U − U ) n ↑ n ↓ = ( U − U ) n n + U − U n ( n −
1) (11)where in the last term we have rewritten the intra-orbitalinteraction for orbital 2 in terms of n = n ↑ + n ↓ insteadof n ↑ and n ↓ . Hence the remaining interaction consistsof an inter-orbital interaction ∼ ( U − U ) and a SSMinteraction ∼ ( U − U ) / U − U is largeror smaller than ( U − U ) /
2, the remaining term is eithera SSM interaction for orbital 2 if U − U < ( U − U ) / U < U ave ≡ ( U + U ) /
2) or an inter-orbital interaction if U − U > ( U − U ) / U > U ave ).We thus identify two new regimes in addition toRegime I ( U ≤ U α ) discussed above: In Regime II the inter-orbital interaction U takes values between thelowest interaction of both intra-orbital interactions andtheir average U ave , i.e. U < U < U ave . This is thecase shown in Fig. 2(a-c). After subtraction of the CIMpotential ∼ U we find that the remaining interaction inRegime II can be written as( U − U ) ˆ n ( ˆ N −
1) + 2( U ave − U ) ˆ n ↑ ˆ n ↓ . (12)Overall this suggests the following decomposition of theHxc potential in Regime II ( U < U < U ave ): v Hxc α [ n ] = v HxcCIM ( U )[ N ]+ v Hxcskew ,α (2( U − U )) [ n ]+ v HxcSSM (2( U ave − U )) [ n ] δ α (13)where δ α is the Kronecker-delta which ensures that theSSM term only contributes to the Hxc potential of orbital2. Note that as U → U ave the SSM term vanishes.On the other hand Regime III occurs when the inter-orbital interaction exceeds the average intra-orbital in-teraction, i.e. U > U ave > U . This was the caseconsidered in Fig. 2(d-f). In this regime the remaininginteraction after subtraction of the CIM ∼ U can berewritten in terms of the Skew interaction (10) and apure inter-orbital interaction part: U − U n ( ˆ N −
1) + ( U − U ave ) ˆ n ˆ n . (14)As can be seen in Fig. 3(d), this inter-orbital term V inter = U ˆ n ˆ n (15)gives rise to a single step at N = 2 which explains theincrease in step height at N = 2 with respect to the CIM,observed in Fig. 2(e,f). Overall this suggests the follow-ing decomposition of the Hxc potential in Regime III ( U > U ave > U ): v Hxc α [ n ] = v HxcCIM ( U )[ N ]+ v Hxcskew ,α ( U − U )[ n ]+ v Hxcinter ( U − U ave )[ n + n ] . (16)We can see that for U = U the Skew part of the Hxcpotential disappears.Hence we have found a decomposition of the Hxc po-tential for a two-orbital model with generic (density-density ) interactions in all three regimes in terms of fourbasic potentials which are shown schematically in Fig. 4.We would like to emphasize at this point that Regime I corresponds to a more natural choice of parameters thanthe other two regimes, since the inter-orbital interaction U is generally smaller than any of the intra-orbital in-teractions U α . Nevertheless, the other regimes might berealized by effective models or possibly by screening ofthe Coulomb interactions. In the next section we willpresent parametrizations of the Hxc potentials in the dif-ferent regimes, making use of its decomposition into thebasic building blocks shown in Fig. 4. B. Generalization of Hxc potential to more thantwo orbitals
For specific choices of parameters we can generalize theHxc potential for the DQD to an arbitrary number oforbitals in a straightforward manner. We concentrate on (a)(b) (c)(d) (e)(f)
FIG. 4. Schematic representation of the four basic Hxcpotentials for building the generic potentials for all threeregimes. (a) Hxc potential for CIM interaction U ˆ N ( ˆ N − Un n .(c,d) Hxc potential for intra-orbital (i.e. SSM) interactions Un α ↑ n α ↓ . (e,f) Hxc for the Skew interaction U ˆ n ( ˆ N − the physically most relevant Regime I ( U α , U β > U αβ ).If we choose the inter-orbital interaction to be constant, U αβ ≡ U (cid:48) , which thus has to be smaller than all of theintra-orbital interactions, U (cid:48) < U α , we can rewrite theinteraction in a similar manner as in Eq. (8) in terms ofa CIM term ∼ U (cid:48) for all the electrons N = (cid:80) α n α andSSM terms ∼ δU α ≡ U α − U (cid:48) for the individual orbitalsas V int = 12 U (cid:48) ˆ N ( ˆ N −
1) + (cid:88) α δU α ˆ n α ↑ ˆ n α ↓ (17)where δU α = U α − U (cid:48) . This suggests to decompose theXC functionals in complete analogy to the two-orbitalcase in Regime I as v Hxc α [ n ] = v HxcCIM ( U (cid:48) )[ N ] + v HxcSSM ( δU α )[ n α ] . (18)In Sec. V B we will see that this decomposition of theHxc potential leads to excellent results for multi-orbitalQDs. For a more general choice of interaction parame-ters, the above decomposition is likely to become morecomplicated. This will be the focus of future work. C. The effect of Hund’s rule coupling
So far we have neglected the effect of Hund’s rule cou-pling on the Hxc potentials. In Fig. 5 we show the sta-bility diagram and the corresponding reverse-engineeredHxc potential for the case of a CIM type direct interac-tion part ( U = U = U ) plus the full Hund’s couplingcontribution ( J H ). Both the stability diagram and thereverse-engineered Hxc potential shown in Fig. 5 resem-ble the ones for the case with U < U α without Hund’scoupling (cf. Fig. 1(d-f)). Only the size of the vertexregions changes in the stability diagram, and correspond-ingly in the Hxc potentials only the step heights change.Moreover, by switching off the spin-flip term in (1) wefind that it does not have any effect on the densities andconsequently on the Hxc potentials and thus can be ne-glected. Hence in the following considerations we onlyneed to take into account the density-density part of theHund’s coupling in (next to last term in Eq. 1).In the spirit of the previous section we can rewrite thedensity-density part of the Hund’s rule coupling term interms of a (negative) CIM interaction and (positive) SSMinteractions for the remaining orbitals plus a remaining positive interaction part: V H = − J H (cid:88) σ ˆ n σ ˆ n σ = − J H ˆ n ˆ n + J H (cid:88) σ ˆ n σ ˆ n σ = − J H N ( ˆ N −
1) + J H (cid:88) α ˆ n α ↑ ˆ n α ↓ + J H (cid:88) σ ˆ n σ ˆ n σ (19)where in the last term ¯ σ denotes the oposite spin of σ .The last term gives rise to a step at N = 2 of height J H in the Hxc potential similar to the inter-orbital inter-action term but with step height J H instead of 2 U (cf.Fig. 3(d)).When adding the density-density contribution of theHund’s rule coupling to the direct interaction part inRegime I ( U ≤ U , U ), we can rewrite the interactionin terms of a CIM interaction, SSM terms, and the lastterm of the Hund density-density interaction (19) as V int = U − J H N ( ˆ N −
1) + (cid:88) α ( δU α + J H )ˆ n α ↑ ˆ n α ↓ + J H (cid:88) σ ˆ n σ ˆ n σ (20)where as defined in the previous section δU α = U α − U . Hence all terms can be modeled by the basic Hxcpotentials shown in Fig. 4: v Hxc α [ n ] = v HxcCIM ( U − J H )[ N ]+ v HxcSSM ( δU α + J H )[ n ]+ v Hxcinter ( J H / n + n ] . (21) D. Parametrization of the basic Hxc potentials
In the zero temperature limit, the steps in the Hxcpotential become infinitely sharp and thus can be de-scribed by simple step functions. At finite temperature,however, the steps are broadened in a non-trivial way.For the SSM at finite temperature an exact expressionfor the Hxc potential can be derived: v HxcSSM [ n ] = U + 1 β ln (cid:32) x + (cid:112) x + e − βU (1 − x )1 + x (cid:33) (22)where x = n − -6-4-2 0 2 -6 -4 -2 0 2 v v (0,0)(1,0)(2,0) (0,1)(0,2)(1,1)(2,1) (1,2)(2,2)(a) 012 0 1 2 n n H H H (b) FIG. 5. Effect of Hund’s rule coupling on (a) Stability dia-gram and (b) Hxc potential of orbital 1 for CIM interactionplus Hund’s rule coupling, V int = U ˆ N ( ˆ N −
1) + V Hund for U = 2 J H . Here due to symmetry the Hxc potential for or-bital 2 can simply be obtained by reflection along the n = n line. All energies in units of J H . other three basic Hxc potentials into which the genericHxc potential can be decomposed, namely the CIM, theInter-orbital, and the Skew potential (see Fig. 4). Westart with the CIM potential and show that an excel-lent parametrization of the Hxc potential can be achievedby simply summing the (exact) SSM potential (22) overthe charging states of the dot, and shifting and rescalingit such that the potential does not become negative orlarger than (2 M − U : v HxcCIM [ N ] = (2 M − Uv maxCIM × M − (cid:88) J =1 (cid:2) v HxcSSM [ N − J + 1] − v HxcSSM [ − J + 1] (cid:3) (23)where v maxCIM = M − (cid:88) J =1 (cid:2) v HxcSSM [2 M − J + 1] − v HxcSSM [ − J + 1] (cid:3) (24)is the maximal value that the sum in (23) aquires at N =2 M . The prefactor (2 M − U/v maxCIM thus rescales thepotential such that the potential yields the exact value(2 M − U at N = 2 M . As can be seen in Fig. 6(a), theagreement with the exact result is quite remarkable.For the inter-orbital potential we find a goodparametrization describing the step at total N = 2 againin terms of the SSM potential, as v Hxcinter ( U, β )[ N ] = v HxcSSM (2 U, β ∗ )[ N/
2] (25)where we have replaced the actual inverse temperature β by an effective reduced value, β ∗ = 0 . β and the stepheight is increased by a factor of 2 compared to the SSM.The agreement with the exact potential is very good ascan be seen in Fig. 6(b).Finally, for the Skew interaction, we parametrize theHxc potential in a similar way as the Hxc potential forthe CIM, by summing two SSM potentials, one for each ofthe steps, and shifting and rescaling so that the potential v H x c parametrizationexact β = 5 β = 10 v H x c a) CIM b) INTERc) SKEW d) SKEW FIG. 6. Comparison of parametrized and exact Hxc poten-tials as a function of N = n + n for three basic interactions:(a) CIM interaction ( U = U = U > U > U = U = 0); (c) Skew interaction( U = 2 U > U = 0); All energies in units of thesmallest non-zero interaction ( U ). does not become negative or larger than the maximumvalue: v Hxcskew ,α ( U )[ n ] = α Uv maxskew (cid:88) J =0 , (cid:8) v HxcSSM ( U ) (cid:2) n + n − J (cid:3) − v HxcSSM ( U )[ − J ] (cid:9) (26)where v maxskew = (cid:88) J =0 , (cid:8) v HxcSSM ( U )[3 − J ] − v HxcSSM ( U )[ − J ] (cid:9) . (27)Also here the agreement with the exact potential is verygood as can be seen in Figs. 6(c+d).We have thus found parametrizations of the four ba-sic Hxc potentials. It should be noted, however, that athigher temperatures the exact CIM and Inter-orbital po-tentials (which in the zero temperature limit only dependon total N ) acquire also a dependence on the difference δN ≡ n − n which has not been taken into accounthere. This discrepancy of our parametrizations might beresponsible for some of the moderate deviations of ourDFT results from the exact ones discussed in the nextsection. V. DFT CALCULATIONS
We are now going to apply the above developedparametrization of the Hxc potential in actual DFT cal-culations. To this end we solve the KS equations whichfor the multi-orbital QD in the GCE are given by: n α = 2 f ( v α + v Hxc [ n ]) for α = 1 . . . M . (28) Since the sharp step features in the Hxc potentials areexpected to prevent the convergence of the usual self-consistency procedure in the density, here we make useinstead of a multidimensional generalization of the bisec-tion approach as proposed before in Ref. 54 for findingthe root of the multidimensional function F α [ n ] ≡ n α − f ( v α + v Hxc [ n ]) . (29)In the following we study the evolution of the density n of multi-orbital QDs as a function of the applied gate v g for different parameter sets. The gate v g exerts a totalshift of the QD levels (cid:15) α and hence the total gate fororbital α is given by v α = (cid:15) α + v g . (30)Consequently, the differences in the gate potentials be-tween different orbitals remain constant as the gate v g changes, δv αβ ≡ v α − v β = (cid:15) α − (cid:15) β . In the following wewill usually take the particle-hole symmetric (phs) pointgiven by (cid:15) ∗ α = − U α − (cid:80) β (cid:54) = α U αβ as the reference system. A. Results for the double quantum dot
We now study the DQD, and start by consideringthe degenerate case in Regime I, i.e. U = U > U where δN = 0. Fig. 7 compares the exact densities withthe ones computed in DFT using the Hxc potential forRegime I, Eq. (9), together with the parametrizations ofthe SSM, Eq. (22), and the CIM, Eq. (23), respectively.We see that the DFT results correctly describe all thefeatures of the densities as a function of gate. At low tem-peratures, the width of the central step (around v g = 0) isgiven by U i while the other two step widths correspond to U . At higher temperatures our parametrization leadsto moderate discrepancies in the slopes of the central stepthat disappear as the temperature approaches zero.Next we consider the situation where the intra-orbitalCoulomb repulsions are different, U > U > U . InFig. 8(a,b), the occupations n i are presented as a func-tion of the gate v g for two different temperatures. Atlow temperatures [Fig. 8(a)] at large negative gate volt-age ( v g < − .
5) both orbitals of the DQD are completelyfilled ( n α ∼ U ) becomes half-filledaround v g ∼ − .
5, and then around v g ∼ − . U ) becomes half-filled. Upon further increase of the gate, the sequence ofemptying is reversed, as first the orbital with the higherinteraction and thus lower gate ( v ) is emptied around v g ∼ . v ) is emptied. At higher temperaturesextra steps develop in the evolution of the density versusgate voltage, as can be seen in Fig. 8(b). The appear-ance of new steps can be understood by the path taken inthe n − n plane as the gate voltage changes, shown inthe inset of Fig. 8(b) for different temperatures. At low0 -2 0 2v g n α β=3β=10β=50 -2 0 2v g n α DFTexact
FIG. 7. Density n = ( n , n ) as function of the gate voltage v g for different temperatures when U = U = 3 U > U . temperatures the path essentially follows three straightline segments, along the lower border, across the planeand finally along the upper border, thus avoiding extrasteps of the CIM potential at N = 1 and N = 3. Asthe temperature increases the path becomes smoother,and passes through the N = 1 and N = 3 steps of theCIM potential, leading to the extra steps in the evolu-tion of the densities at higher temperature. While forlow temperatures the agreement of the DFT results withthe exact ones is excellent, at higher temperatures de-viations appear. Although DFT qualitatively capturesthe appearance of the extra steps in the evolution of thedensity versus gate voltage, their heights are not cor-rectly reproduced in DFT. Presumably this discrepancycan be attributed to the development of a δN -dependenceof the CIM potential at finite temperature, and will beaddressed in future work.Finally, we turn our attention to Regimes II and III,which are both characterized by the appearance of thepeculiar “Skew” term in the Hxc potential. Fig. 9(a) di-rectly compares the evolution of the density as a functionof the gate in both regimes. As we can see the behaviouris actually quite similar for both regimes, and not so dif-ferent from Regime I (cf. Fig. 8): As the gate increases,first the orbital with the higher interaction (here U ) be-comes half-filled, and then the orbital with the lower in-teraction ( U ). Then upon further increase of the gate,the order of emptying is reversed. Due to the higherinter-orbital interaction in Regime III, the width of thecentral plateau is increased for both orbitals. Again, atlow temperature the agreement with the exact results isexcellent, but at higher temperatures moderate quanti-tative deviations occur (not shown).In order to investigate the influence of the Skew termin the Hxc potential on the evolution of the densities,we next concentrate on Regime II and explore differentpaths in the n − n plane. To this end we fix the energysplitting δv = v − v between the orbitals to different -2 0 2v g n n -2 0 2v g n α DFTexact n n β =1 β =6 β =10 β =50 a) b) β=6β=50 FIG. 8. Density n = ( n , n ) as a function of the gate volt-age v g for U = 3 U , U = 2 . U (Regime I) for (a) lowand (b) high temperatures. The inset of panel (b) shows thepath in the n − n plane as the gate is varied for differenttemperatures. The grey lines show the steps of the CIM andSSM terms that appear in the Hxc potentials. All energies inunits of the smallest interaction U . -5 0 5v g n α n RIIn RIIn RIIIn RIII -10 -5 0v g n δ v=0n δ v=0n δ v=-1n δ v=-1 n n DFT -10 -5 0v g DFTexact a) b)
FIG. 9. (a) Comparsion of the evolution of the density n =( n , n ) with the gate voltage v g in Regime II ( U = 2 U =4 U ) and Regime III ( U = U = 4 U ). (b) Comparisonof density evolution for different two different values of thesplitting δv in Regime II ( U = 2 U ). The inset shows thedifferent paths in the density plane. β = 20 /U everywhere.All energies in units of U . values while the total gate changes, i.e. v = δv + v g and v = v g . Fig. 9(b) shows the evolution of the density fortwo different values of δv and correspondingly differentpaths in the n − n plane (shown in the inset). For δv = 0we observe an interesting effect. As the gate increases,the occupation of orbital 2 decrease in two steps, firstto half filled and then further to zero, while the firstorbital remains fully occupied. Then around v g = − n α -5 0 5v g n α DFTexact -10 -5 0v g n n n n a) b)c) d) FIG. 10. Local occupations for the triple (a, b) and quadruple(c, d) QD as function of the gate voltage. In the left panels theQD levels are taken at particle-hole ( v α = (cid:15) ∗ α + v g ) and in theright the impurities level is set to zero ( v α = v g ). β = 20 /U (cid:48) and U = 5 U (cid:48) , U = 4 U (cid:48) , U = 3 U (cid:48) , U = 2 U (cid:48) . All energiesin units of U (cid:48) . again to quarter filling, n = n = ∼ .
5. This non-monotonic behaviour of the occupation of orbital 2 isreminiscent of the so-called level occupation switching(LOS) . We find similar behaviour in Regime III (notshown).
B. Results for more than two orbitals
Finally we apply the generalization of the Hxc poten-tial (18) for more than two orbitals to DFT calculationsof multi-orbital QDs. Fig. 10 shows the evolution of thedensity n as a function of the applied gate voltage v g forthree (a,b) and four-level (c,d) QD with all intra-orbitalCoulomb repulsions U α different and constant interdotrepulsion U (cid:48) ( U > U > . . . > U (cid:48) ) at low temperature.In panels (a) and (c) the gate v g is applied w.r.t. thephs point, i.e. (cid:15) ∗ α = − U α − (cid:80) β (cid:54) = α U αβ . In this casethe path in the three- or four-dimensional density spaceavoids the steps in the CIM Hxc potential away fromhalf-filling ( N = M ) resulting in only three plateaus inthe density evolution with the gate, in a similar way as inthe DQD [cf. Fig. 8(a)]. On the other hand, in panels (b)and (d) (where (cid:15) α = 0 and thus δv αβ = 0) two (three) ex-tra steps related to the inter-orbital Coulomb repulsionsappear in the triple (quadruple) QD. The agreement be-tween the DFT and the exact results is remarkable in allcases, showing that the generalization of Eq. (18) of theHxc potential to more than two orbitals is valid. Findingsimilar expressions for a more general choice of parame-ters will be the focus of future work. VI. CONCLUSIONS
In this work we have obtained Hxc potentials for dou-ble quantum dots in the grand-canonical ensemble sub- ject to generic density-density interactions and Hund’srule coupling by reverse-engineering from exact many-body solutions. The structure of the Hxc potentials con-sists of ubiquitous steps whose exact positions depend onthe regime defined by the interaction parameters. Thisstructure can be understood and derived from an anal-ysis of the stability diagrams. In a second step we wereable to rationalize the step structure of the Hxc potentialby a decomposition of the interaction into basic compo-nents. This decomposition allows to write the Hxc po-tential of the system as a sum over basic Hxc poten-tials, which can be parametrized in a straightforwardmanner. Importantly, the decomposition into basic po-tentials can be generalized to multi-orbital systems withmore than two orbitals. DFT calculations employing thethus parametrized Hxc potentials for double, triple andquadruple quantum dots show excellent agreement withexact results at low temperatures. At higher temper-atures, we find moderate quantitative deviations fromthe exact results that we attribute to the modulation ofstep widths for finite δn = n − n not captured in ourparametrization.The parametrization of the Hxc potential derived herecould be directly applied e.g. to the description of singleatoms where density-density and Hund’s coupling are thedominating terms of the Coulomb interaction. Possiblefurther applications regard the description of transportthrough multi-orbital quantum dots or molecules coupledto leads. Due to the similarity between broadening byfinite temperature on the one hand and finite couplingto the leads on the other hand, we expect that the Hxcpotentials for finite coupling to the leads have a simi-lar structure to the ones discussed here. One way toincorporate finite coupling to the leads in the Hxc po-tential is by introduction of an effective temperature. For the description of non-equilibrium effects the i-DFTframework may be employed which in addition to theHxc gate potential requires a parametrization of the xcbias. This would also allow one to compute many-bodyspectral functions of interacting multi-orbital systems. ACKNOWLEDGMENTS
We acknowledge funding by the grant Grupos Con-solidados UPV/EHU del Gobierno Vasco (IT1249-19) aswell as the grant of the Ministerio de Econom´ıa, Indus-tria y Competitividad, Gobierno de Espa˜na (MINECO) -Agencia Estatal de Investigaci´on (FIS2016-79464-P) andEuropean Regional Development Fund (FEDER), Euro-pean Union.2 ∗ [email protected] P. Hohenberg and W. Kohn, Phys. Rev , B864 (1964). W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965). R. M. Dreizler and E. K. U. Gross,
Density FunctionalTheory (Springer, Berlin, 1990). J.P. Perdew, Phys. Rev. Lett. , 1665 (1985). A.D. Becke, Phys. Rev. A , 3098 (1988). J.P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. , 3865 (1996); ibid. , 1396 (1997)(E). J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov,G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke,Phys. Rev. Lett. , 136406 (2008). A.D. Becke, J. Chem. Phys. , 5648 (1993). J. Heyd, G. Scuseria, and M. Ernzerhof, J. Chem. Phys. , 8207 (2003). J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Phys.Rev. Lett. , 1691 (1982). E. Sagvolden and J. P. Perdew, Phys. Rev. A , 012517(2008). P. Gori-Giorgi and A. Savin, Int. J. Quantum Chem. ,2410 (2009). L. J. Sham and M. Schl¨uter, Phys. Rev. B , 3883 (1985). P. W. Anderson, Phys. Rev. , 41 (1961). J. Hubbard, Proc. R. Soc. Lond. A , 238 (1963). A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,Rev. Mod. Phys. , 13 (1996). G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O.Parcollet, and C. A. Marianetti, Rev. Mod. Phys. , 865(2006). D. Jacob, K. Haule, and G. Kotliar, Phys. Rev. B ,195115 (2010). C. Weber, D. J. Cole, D. D. ORegan, and M. C. Payne,Proc. Nat. Acad. Sci. , 5790 (2014). M. Karolak, G. Ulm, T. O. Wehling, V. Mazurenko, A.Poteryaev, and A. Lichtenstein, J. Electron Spectrosc. Re-lat. Phenom. , 11 (2010). R. Requist and E. K. U. Gross, Phys. Rev. B , 125114(2019). J. P. Coe, Phys. Rev. B , 165118 (2019). L. Mazouin, M. Sauban`ere, and E. Fromager, Phys. Rev.B , 195104 (2019). K. Haule, Phys. Rev. Lett. , 196403 (2015). O. Gunnarsson and K. Sch¨onhammer, Phys. Rev. Lett. ,1968 (1986). K. Sch¨onhammer, O. Gunnarsson, and R. M. Noack, Phys.Rev. B , 2504 (1995). N. A. Lima, L. N. Oliveira, and K. Capelle, Europhys.Lett. , 601 (2002). N. A. Lima, M. F. Silva, L. N. Oliveira, and K. Capelle,Phys. Rev. Lett. , 146402 (2003). R. L´opez-Sandoval and G. M. Pastor, Phys. Rev. B ,035115 (2003). G. Xianlong, M. Polini, B. Tanatar, and M. P. Tosi, Phys.Rev. B , 161103(R) (2006). K. Capelle and V. L. Campo Jr., Phys. Rep. , 91 (2013). V. Brosco, Z.-J. Ying, and J. Lorenzana, Sci. Rep. , 2172(2013). Z.-J. Ying, V. Brosco, and J. Lorenzana, Phys. Rev. B ,205130 (2014). D. Carrascal, J. Ferrer, J. Smith, and K. Burke, J. Phys.:Condens. Matter , 019501 (2015). T. M¨uller, W. T¨ows, and G. M. Pastor, Computation ,66 (2019). C. Verdozzi, G. Stefanucci, and C.-O. Almbladh, Phys.Rev. Lett. , 046603 (2006). C. Verdozzi, Phys. Rev. Lett. , 166401 (2008). S. Kurth, G. Stefanucci, E. Khosravi, C. Verdozzi, andE. K. U. Gross, Phys. Rev. Lett. , 236801 (2010). J. I. Fuks and N. T. Maitra, Phys. Rev. A , 062502(2014). J. I. Fuks and N. T. Maitra, Phys. Chem. Chem. Phys. ,14504 (2014). N. Dittmann, J. Splettstoesser, and N. Helbig, Phys. Rev.Lett. , 157701 (2018). N. Dittmann, N. Helbig, and D. M. Kennes, Phys. Rev. B , 075417 (2019). G. Stefanucci and S. Kurth, Phys. Rev. Lett. , 216401(2011). J. P. Bergfield, Z.-F. Liu, K. Burke, and C. A. Stafford,Phys. Rev. Lett. , 066801 (2012). P. Tr¨oster, P. Schmitteckert, and F. Evers, Phys. Rev. B , 115409 (2012). M. Seidl, J. P. Perdew, and M. Levy, Phys. Rev. A , 51(1999). A. Mirtschink, M. Seidl, and P. Gori-Giorgi, Phys. Rev.Lett. , 126402 (2013). A. Georges, L. de’ Medici and J. Mravlje, Ann. Rev. Con-dens. Matt. Phys. , 137 (2012). N. Mermin, Phys. Rev. , A1441 (1965). G. Stefanucci and S. Kurth, Phys. Stat. Sol. (b) , 2378(2013). S. Kurth and G. Stefanucci, J. Phys.: Condens. Matter ,413002 (2017). T. Dimitrov, H. Appel, J. Fuks, and A. Rubio, New J.Phys. , 083004 (2016). E. Perfetto and G. Stefanucci, Phys. Rev. B , 081409(R)(2012). G. Xianlong, A.-H. Chen, I. V. Tokatly, and S. Kurth,Phys. Rev. B , 235139 (2012). P. G. Silvestrov and Y. Imry, New J. Phys. , 125 (2007). Y. Kleeorin and Y. Meir, Phys. Rev. B , 045118 (2017). S. Kurth and G. Stefanucci, Phys. Rev. B , 241103(R)(2016). N. Sobrino, R. D’Agosta, and S. Kurth, Phys. Rev. B ,195142 (2019). G. Stefanucci and S. Kurth, Nano Lett. , 8020 (2015). D. Jacob and S. Kurth, Nano Lett.18