Table Of ContentCorrelation energy of two-dimensional systems:
Toward non-empirical and universal modeling
S. Pittalis,1,2,∗ E. Ra¨s¨anen,1,2,3,† C. R. Proetto,1,2,‡ and E. K. U. Gross1,2
1Institut fu¨r Theoretische Physik, Freie Universit¨at Berlin, Arnimallee 14, D-14195 Berlin, Germany
2European Theoretical Spectroscopy Facility (ETSF)
3Nanoscience Center, Department of Physics, University of Jyv¨askyl¨a, FI-40014 Jyv¨askyl¨a, Finland
(Dated: January 30, 2009)
The capability of density-functional theory to deal with the ground-state of strongly correlated
low-dimensionalsystems,suchassemiconductorquantumdots,dependsontheaccuracyoffunction-
9
alsdevelopedfortheexchangeandcorrelation energies. Hereweextendasuccessfulapproximation
0
for the correlation energy of the three dimensional inhomogeneous electron gas, originally intro-
0
duced by Becke [J. Chem. Phys. 88, 1053 (1988)], to the two-dimensional case. The approach
2
aimstonon-empiricalmodelingofthecorrelation-holefunctionssatisfyingasetofexactproperties.
n Furthermore, the electron current and spin are explicitly taken into account. As a result, good
a performance is obtained in comparison with numerically exact data for quantumdotswith varying
J externalmagnetic field, and for the homogeneous two-dimensional electron gas, respectively.
0
3 PACSnumbers: 73.21.La,71.15.Mb
]
el I. INTRODUCTION plies onlyto closed-shellsystems withzerospinandzero
- current. Inadditiontothecorrelation-energyfunctional,
r
t newexchange-energyfunctionalshavebeendevelopedfor
s Density-functional theory1 (DFT) maps the compli- finite 2D systems in our foregoing works.9–11
.
t cated many-particle problem onto a simple one of non- In this paper we develop a correlation-energy func-
a
interacting electrons moving in an effective, local poten-
m tional in 2D. In the derivation, along the lines of the
tial, the Kohn-Sham (KS) potential. The latter is con- work of Becke12 for the 3D case, we introduce a model
-
structed in such a way that the ground-state density of
d forspin-dependentcorrelation-holefunctionssatisfyinga
n thenon-interactingparticlesreproducestheground-state setof exactproperties in 2D.As a result, we find a spin-
o density of the interacting system. Practical success of
andcurrent-dependentapproximationforthecorrelation
c the approach depends on finding good approximations
energy. Incomparisonwithnumericallyexactresults,the
[
for the exchange-correlation (xc) energy which, through obtained accuracy is found to be superior to the LSDA.
2 functional derivation with respect to the particle den- However,wealsofind,andelucidate,thatfurthermodel-
v sity, defines the xc part of the KS potential. Most of the
ing of the dependency on the averageelectron density of
3 approximations developed so far have focused on three-
those parameters describing the size of the correlation-
8 dimensional (3D) systems, i.e., atoms, molecules, and
hole functions in terms of the size of the exchange-hole
2
solids. Many advanceshave been made beyondthe com-
4 (x-hole) functions would be required. The applications
monly used local(spin) density approximation[L(S)DA]
. of the functional to a set of few-electron QDs with var-
0 by means of, e.g., generalized gradient approximations,
ious relative amounts of correlation, ground-state spins,
1 orbital functionals, and hybrid functionals.2 Such efforts
8 electron currents, and external magnetic fields, confirm
for two-dimensional (2D) systems have been relatively
0 the overallusefulness of the approach.
scarce despite the rapidly increasing experimental and
:
v theoretical interest in quasi-2D structures such as semi-
i
X conductor layers and surfaces, quantum-Hall systems, II. THEORY
graphene, and various types of quantum dots3 (QDs).
r
a When using DFT, QDs are most commonly treated Within the Kohn-Sham (KS) method of spin-DFT,13
using the 2D-LSDA exchange4 combined with the 2D- the ground-state energy and spin densities ρ↑(r) and
LSDA correlation parametrized first by Tanatar and ρ↓(r) of a system of N = N↑+N↓ interacting electrons
Ceperley5 and later, with more satisfactory spin depen- are determined. The total energy of the interacting sys-
dence, by Attaccalite et al.6 In many cases, the LSDA tem is written as a functional of the spin densities14
(prefix “2D” omitted below) performs relatively well
compared,e.g.,withquantumMonteCarlocalculations.7 Ev[ρ↑,ρ↓] = Ts[ρ↑,ρ↓]+ drv(r)ρ(r)
Nevertheless,thereisalackof2Dfunctionalstodealwith Z
+ E [ρ]+E [ρ ,ρ ] (1)
diversefew-electronQDsystems,especiallyinthestrong- H xc ↑ ↓
correlationregime. Onlyveryrecently,alocalcorrelation where T [ρ ,ρ ] is the kinetic energy functional of non-
s ↑ ↓
functional was developedin 2D within the Colle-Salvetti interacting electrons with spin densities ρ , ρ . v is
↑ ↓
approach, which was found to outperform the LSDA.8 (at vanishing external magnetic field) the external (lo-
In its current form, however, this local functional ap- cal) scalar potential acting upon the interacting system,
2
E [ρ] is the classical electrostatic or Hartree energy of sum rule
H
•
the total charge density ρ = ρ +ρ , and E [ρ ,ρ ] is
↑ ↓ xc ↑ ↓
thexcenergyfunctional. Excmaybefurtherdecomposed dss¯hσc,σλ′(r,s)=0; (8)
into the exchangeenergy,E , andcorrelationenergyE .
x c Z
We have already considered E in Ref. 9, and thus we
x
here focus on E . Our starting point is the formal ex- correct short-range behavior for s 0;
c • →
pression for E in terms of the correlation-hole (c-hole)
c a proper decay in the limit s ;
function, • →∞
characteristic size assumed to be proportional to
Ec[ρ↑,ρ↓]= 21 dr1 dr2 rρσ(r1r) hσcσ′(r1,r2), • the characteristicsizeofthecorrespondingangular
σσ′ Z Z | 1− 2| averagedx-hole.
X
(2)
where The short-range behavior of the c-hole can be worked
out by considering the electronic cusp conditions for the
1
hσσ′(r ,r )= dλhσσ′(r ,r ) 2D electronic wave function.15 The angular average of
c 1 2 c,λ 1 2
Z0 the 2BRDM in Eq. (6) is then given by
1
=Z0 dλhσλσ′(r1,r2)−hσx(r1,r2)δσσ′ (3) P¯2σ,σλ(r,s→0)≈Aσσ(r)s2 1+ 32λs , (9)
is the c-hole function. Here σσ′ = , , , (cid:18) (cid:19)
{ } {↑↑ ↓↓ ↑↓ ↓↑}
and the parameter λ [0,1] is the electronic coupling and
strength.1Intheabove∈expression,thec-holecorresponds
to the full hole function hσλσ′ subtracted by the x-hole P¯2σ,σλ¯(r,s→0)≈Aσσ¯(r)(1+2λs), (10)
function defined as
for the same-spin and opposite-spin elements, respec-
hσ(r ,r )= | Nk=σ1ψk∗,σ(r1)ψk,σ(r2)|2. (4) tively. Notethatthesymbolσ¯ alwaysdenotesspinoppo-
x 1 2 − ρ (r ) site to σ, whereas σ′, when used, can be equal to either
P σ 1
Note that here hσ is defined to have a negative sign in σ or σ¯. Here the coefficients Aσσ′(r) have a spatial de-
x pendence. Although the cusp condition as developed in
contrast to the standard definition [see, e.g., Eq. (4) in
Ref. 15 pertains to the center-of-mass and relative coor-
Ref. 9]. The λ-dependent hole function is given by dinate, it can be shown, as noted by Becke,12 that, at
Pσσ′(r ,r ) the order considered above, the same results apply also
hσλσ′(r1,r2)= 2,ρλσ(r11) 2 −ρσ′(r2), (5) forTthheeschooorrtd-riannagteebsyeshtaevmioorfofrtahned2Ds.x-holeisknown.9,16
where Pσσ′ is the two-body reduced density matrix Thus, it is possible to consider the following model for
2,λ thesame-spinandopposite-spinc-holefunctions,respec-
(2BRDM) 12
tively:
P2σ,σλ′(r1,r2) = N(N −1) d3 d4... dN h¯σσ(r,s) = B (r) D (r)+ 2λB (r)s s2
× Ψ∗λ(1,2,...,ZN)ΨZλ(1,2,.Z..,N). (6) c,λ F(cid:20) (σγσ (r)−s),σ 3 σσ (cid:21) (11)
σσ
Here, Ψ (1,2,...,N) stands for the ground-state many- ×
λ
bodywavefunctionwhichistheexactsolutionoftheelec- and
tronicsystemwithaCoulombcouplingstrengthλ, dN
denotes the spatial integrationand spin summation over h¯σσ¯(r,s) = [B (r) ρ (r)+2λB (r)s]
R c,λ σσ¯ − σ¯ σσ¯
the N-th spatial spin coordinates (rN,σN), and we have F(γσσ¯(r)s). (12)
identified σ = σ, σ = σ′. Note that the spin densities ×
1 2
in Eq.(5) arethe same as in the actual,fully interacting HereB :=A /ρ andB :=A /ρ arecoefficients
σσ σσ σ σσ¯ σσ¯ σ
system. to be determined, and
As it is well known, in determining the correlationen-
ergy it is sufficient to know the angular average of the 1 1( ρ )2 j2
c-hole function. In 2D is natural to consider the cylin- Dσ := τσ ∇ σ p,σ ; (13)
2 − 4 ρσ − ρσ !
drical average,given by
h¯σσ′(r,s)= 1 2πdφ hσσ′(r,r+s), (7) where τσ = Nk=σ1|∇ψk,σ|2 is the (double
c,λ 2π Z0 c,λ of the) kinetic-enPergy density, and jp,σ =
where r1 = r, r2 = r+s, and φ is the angle between r 21i Nk=σ1 ψk∗,σ(∇ψk,σ)− ∇ψk∗,σ ψk,σ is the spin-
and s. The modeling of the c-hole functions is based on dependenht paramagnetic(cid:16)current(cid:17)densiity. In Eqs. (11)
P
satisfying the following conditions: and (12), the functions F(γσσ′(r)s) are introduced to
3
ensure the decay of the c-holes in the limit s . We and
choose them to have the form F(x) = exp( x→2),∞which
− ρ (r)[s z (r)]F (γ (r)s)
seems appropriate in the case of finite 2D systems. The ¯hσσ¯(r,s) = σ¯ − σσ¯ σσ¯
parameters γσσ′(r) are determined by the zero-integral c 2zσ2σ¯(r)
constraint of Eq. (8) (see below). [2z (r) ln(2z (r)+1)] . (21)
σσ¯ σσ¯
× −
Next, we introduce the characteristic sizes of the c-
holes, z and z , for which the corresponding c-hole Finally, we enforce the sum rules in Eq. (8) giving
σσ σσ¯
function vanishes, i.e.,
3√π
γ (r)= , (22)
h¯σc,σλ(r,zσσ)=0; h¯σc,σ¯λ(r,zσσ¯)=0. (14) σσ 4zσσ(r)
Implying these conditions, the coefficients B and B and
σσ σσ¯
in Eqs. (11) and (12) can be written as
√π
γ (r)= . (23)
D ρ σσ¯ 2z (r)
B = σ , B = σ¯ . (15) σσ¯
σσ 1+ 32λzσσ σσ¯ 1+2λzσσ¯ This concludesthe derivationofthe c-holefunctions h¯σσ
Itshouldbenoted,however,thatthesizeparameterszσσ and ¯hσσ¯, apart from the determination of constants cσσ,
and zσσ¯ are functions of r. As for the 3D functional, we and cσσ¯, respectively.
assume the size of the c-hole to be proportional to the From the c-hole functions we can calculate the c-hole
size of the x-hole. Thus, we set potentials as
∞
zσσ(r) := cσσ |Uxσ(r)|−1+|Uxσ(r)|−1 Ucσσ′(r)=2π dshσcσ′(r,s). (24)
= 2cσσ(cid:2)|Uxσ(r)|−1, (cid:3) (16) Z0
For the same- and opposite-spincases ofour approxima-
and
tion, we find respectively
zσσ¯(r):=cσσ¯ |Uxσ(r)|−1+|Uxσ¯(r)|−1 , (17) Uσσ(r) = 16 (8 3π)D (r)z2 (r)
where Uσ is the x-hole(cid:2)potential9,17 for spin (cid:3)σ, and c c 81π − σ σσ
x σσ 2
and c are constants to be determined (see below). As 2z (r) 3ln z (r)+1 , (25)
σσ¯ σσ σσ
suggestedbyBecke,12theproportionalityisplausibledue × (cid:20) − (cid:18)3 (cid:19)(cid:21)
to the fact that each electron is surrounded by its Fermi
and
hole,andhencetheelectrostaticinteractionbetweentwo
electrons can be expected to be screened beyond some Uσσ¯(r) = (2 π)ρ (r)
characteristic length proportional to the average size of c − σ¯
[2z (r) ln(2z (r)+1)]. (26)
thex-hole. Theargumentisnotspin-related,beyondthe × σσ¯ − σσ¯
fact that the characteristic length for the x-hole could
The correlation energies are given by
be different for spin up than for spin down. Of course,
this is only a simplification, but the physical picture is Eσσ′ = 1 drρ (r)Uσσ′(r). (27)
appealing, and it has led to good results in 3D atomic c 2 σ c
systems. Z
The λ-dependent c-hole functions can now be written Thus
as
E [ρ ,ρ ]=E↑↑+E↓↓+2E↑↓, (28)
c ↑ ↓ c c c
2 s z (r)
h¯σσ(r,s)= λD (r) − σσ s2F(γ (r)s),
c,λ 3 σ 1+ 2λz (r) σσ where we have used the condition E↑↓ =E↓↑.
(cid:20) 3 σσ (cid:21) c c
(18) Alternatively, we can compute the correlation energy
and directly from the c-hole functions. From Eqs. (2) and
(7) we get
s z (r)
h¯σσ¯(r,s)=2λρ (r) − σσ¯ F(γ (r)s). (19)
c,λ σ¯ (cid:20)1+2λzσσ¯(r)(cid:21) σσ¯ Ec[ρ↑,ρ↓]=π dr ds ρσ(r)h¯σcσ′(r,s). (29)
Integrating over λ yields Xσσ′ Z Z
WeremindthatD (r)introducedinEq.(13)vanishes
D (r)[s z (r)]s2F (γ (r)s) σ
¯hσσ(r,s) = σ − σσ σσ for all the single-particle (N = 1) systems.12 There-
c 2z2 (r)
σσ fore, the c-hole and thus Ec vanish as well. In other
2 words, our approximation for the correlation energy is
2z (r) 3ln z (r)+1 , (20)
× σσ − 3 σσ self-interaction free for N =1. 26
(cid:20) (cid:18) (cid:19)(cid:21)
4
TABLE I: Total energies from the full configuration- TABLE II: Similar to Table I but for unpolarized (Sz = 0)
interaction calculations (Ref. 21) for totally spin-polarized quantum dots. The correlation energies from our functional,
(Sz =N/2)quantumdots,exact-exchangetotalenergies,the Ecmodel, have been calculated using the fixed average values
referencecorrelationenergies(Ecref =Etot−EtEoXtX),cσσ yield- cσσ =1.32 and cσσ¯ =0.75.
icnσgσE=cm1o.d3e2l,=anEdcretfh,eELcmSoDdeAlocbotrarienleadtiownithenaerfigxye.daveragevalue N ω Etroetf EtEoXtX Ecref cσσ¯ Ecmodel EcLSDA
2 1 3∗ 3.162 -0.162 0.72 -0.171 -0.199
N ω Etroetf EtEoXtX Ecref cσσ Ecmodel EcLSDA 2 1/4 0.9324† 1.046 -0.114 0.82 -0.102 -0.139
3 1/4 2.081 2.103 -0.0226 1.27 -0.0245 -0.0538 2 1/16 0.3031† 0.373 -0.070 0.96 -0.053 -0.085
3 1/16 0.6908 0.7075 -0.0167 1.41 -0.0144 -0.0382 6 1/4 6.995† 7.391 -0.396 0.73 -0.406 -0.457
6 1/4 7.233 7.296 -0.0640 1.24 -0.0750 -0.1125 12 1/1.892 25.636‡ 26.553 -0.917 0.71 -0.983 -1.000
6 1/16 2.553 2.599 -0.0458 1.36 -0.0441 -0.0795
∗ Analytic solution by Taut from Ref. 18.
† CI data from Ref. 21.
III. NUMERICAL RESULTS ‡ Diffusion QMC data from Ref. 23.
The first task in the numerical applications is to com-
plete the correlation functional by finding approxima-
tions for constants c and c in Eqs. (16) and (17). 98 QMC
σσ σσ¯
EXX
For this purpose, we consider a set of harmonically con-
model
fined QDs, where the external confinement is given by
v(r) = ω2r2/2. Reference results for the correlation en- 96 model*
V) LSDA
ergies can be obtained from e
m
Eref =Eref EEXX, (30) Ω ( 94
c tot− tot − 6
where Eref is the exacttotal energyobtained, e.g., from, ot
tot E t92
an analytic, accurate configuration-interaction (CI) or
quantum Monte Carlo (QMC) calculation, and EXX
refers to the exact exchange. Here we have calculated
90
the EXX energies in the Krieger-Li-Iafrate19 (KLI) ap-
proach19intheoctopusDFTcode.20 Theself-consistent
EXXresult–thex-holepotential,(spin)density,kinetic- 0 2 4 6 8
B (T)
energydensity,andcurrentdensity– isusedasinput for
our correlationfunctional.
FIG. 1: (color online). Total energies (minus the confine-
Table I shows the results for the same-spin case. Now
ment energy) for a six-electron quantumdot as a function of
the QDs are completely spin-polarized with S = N/2.
z themagneticfield. Resultsareshown for thereferencequan-
For each QD we show that value of c which yields the
σσ tumMonteCarlocalculationsfromRef.7(redsolidline),for
referencecorrelationenergy,i.e., Ecmodel =Ecref. Inthese theexact-exchangedata(bluedashedline),forourfunctional
exampleswe find cσσ =1.24...1.41. Thus,the variation with cσσ = 1.32 and cσσ¯ = 0.75 (circles), for our functional
of cσσ is rather smallin view of the factthat the density withcσσ =1.1andcσσ¯ =0.7(crosses),andforthelocal-spin-
parameter,definedinharmonicQDsasr =N−1/6ω−2/3 density approximation (green dotted line). The total energy
s
(Ref. 22), varies from1.9 to 5.3. The secondlastcolumn iscalculatedfromourfunctionalasEtmotodel=EtEoXtX+Ecmodel.
ofTableIshowsthecorrelationenergyobtainedbyusing
a fixed average value of c = 1.32. This leads to the
σσ
maximumdeviationof 23%fromEref. Incomparison, the LSDA correlation is only about 9%. This is in line
∼ c
theself-consistentLSDAcorrelationenergy(lastcolumn) withthewell-knownfactthatboththeL(S)DAexchange
deviates from the reference result by up to 130%. andcorrelation,respectively,becomemoreaccuratewith
Table II shows the results for a set of unpolarized increasing particle number.
(S = 0) QDs in the range 1.5 < r < 5.7. Note that In Fig. 1 we consider the total energy of a more gen-
z s
for N > 2 both same- and opposite-spin components of eral case: A six-electron QD as a function of the mag-
the correlationarepresent,andwehavefixedc =1.32 netic field B directed perpendicular to the dot plane.
σσ
according to the conclusions above. Fixing c = 0.75 Increasing the field leads to non-trivial changes in the
σσ¯
yields deviations of only . 10% from Eref, except for ground-state quantum numbers (S ,L ) and hence to
c z z
the highlycorrelatedcaseofN =2andω =1/16,which “kinks” in the ground-state total energy as a function
showsadeviationof25%. TheLSDAisstillconsiderably of B. As the reference data, we use here the variational
further off the reference result than our functional, but QMC results (red solid line) given in Refs. 7 and 25 for
it performs relatively better than in the polarized case a wide range of B up to total spin polarization. The
discussedabove. Inparticular,whenN =12theerrorin confinement strength is here ω = 0.42168, correspond-
5
ing to a typical confinement of 5 meV when modeling 0
QDs in GaAs.3 Note that the total confinement energy,
i.e., 6Ω=6 ω2+ω2/4, where ω =B/c, has been sub-
c c
tractedfromthe totalenergiestoclarifythecomparison. −0.05
p spin−polarized 2DEG
We point out that the variational QMC method gives
an upper bound for the true total energy. On the ba- u.) −0.1 spin−unpolarized 2DEG
sis of previous comparisons between the variational and a.
diffusion QMC, and exact diagonalization,24 our refer- N (
ence data in Fig. 1 can be expected to overestimate the / c−0.15
E
exact total energy by at most 0.2...0.3 meV. The max-
imum possible errors are smaller in the polarized regime
(B &5) T. −0.2 LSDA
Overall,Fig.1showsreasonableagreementofourfunc- model
tional(circles)withtheQMCdatathroughthefullrange
ofthe magneticfield. However,the functionalyieldssys- −0.25 0 5 10 15 20
tematicallyslightlytoolowcorrelationenergies,andthus rs (a.u.)
too low total energies, even if the possible overestima-
tion of the total energy given by the variationalQMC is FIG. 2: (color online). Correlation energy per electron in a
taken into account. We point out that the functional is homogeneous two-dimensional electron gas of full spin polar-
ization (upper curves) and zero polarization (lower curves).
here applied with the fixed parameters c = 1.32 and
σσ
The solid lines show the result from our functional with the
c = 0.75 suggested by the results in Tables I and II,
σσ¯ original average parametervaluescσσ =1.32 and cσσ¯ =0.75.
which correspond to considerably weaker confining po-
Thedashedlinesshowtothelocalspin-densityapproximation
tentials (smaller values of ω). Obviously, this difference
forthecorrelation,correspondingtothenumericallyexactre-
implies a different average electron density, and thus a sult in this system.6
different rangeof the relativecorrelationenergy. In fact,
if the parameter values are reduced to c = 1.1 and
σσ
c = 0.7, excellent agreement with QMC is found (see
σσ¯ IV. CONCLUSIONS
the crosses in Fig. 1). Hence, it seems that particularly
high precision of our functional would require model-
ing of c and c as a function of the particle density. We have derived a spin- and current-dependent ap-
σσ σσ¯
Most importantly, however, the systematic performance proximation for the correlation energy of finite two-
of our approximation(s) in Fig. 1 demonstrates that the dimensional electron systems. The core of the deriva-
magnetic-fieldeffects,electroncurrents,andspinarecor- tion is a model for the correlation-hole function of both
rectly accounted for. same-spinandopposite spinpairs,respectively,thatsat-
We note thatthe goodaccuracyofthe LSDAin terms isfies a set of exact properties. The excellent results ob-
of total energies (see the dotted line in Fig. 1) is due to tainedforfew-electronquantumdotswithdifferentspin-
thecompensationofrespectiveerrorsintheexchangeand polarization, current, external magnetic field, and cov-
correlation energies.25 On the other hand, we tested our ering a wide range of correlation energies, strongly rec-
correlation functional for the important limit of the ho- ommendfurtherdevelopmentsalongtheconstructionwe
mogeneous2Delectrongas(2DEG),forwhichtheLSDA have presented here.
correlation is exact, and found reasonable agreement as
a function of r for both zero and full spin polarization
s
(seeFig.2). Hereweusedtheoriginalaverageparameter
valuesc =1.32andc =0.75. Notethataccordingto
σσ σσ¯
Acknowledgments
the numerically exact results in Ref. 6, the ground-state
of the 2DEG is unpolarized for 0<r .26.
s
Finally, we point out, that in principle a given func- WethankAriHarjuforthevariationalquantumMonte
tional should be evaluated with KS orbitals obtained Carlo data. This work was supported by the Deutsche
from self-consistent calculation instead of a post-hoc Forschungsgemeinschaft,theEU’sSixthFrameworkPro-
manneraswehavedoneinthiswork. However,thevari- gramme through the Nanoquanta Network of Excellence
ational nature of DFT implies that if one evaluates the (NMP4-CT-2004-500198), and by the Academy of Fin-
total energy with densities which slightly differ from the land. C. R. P. was supported by the European Commu-
self-consistent one, the resulting change in the energy is nity through a Marie Curie IIF (MIF1-CT-2006-040222)
of second order in the deviation of the densities. and CONICET of Argentina through PIP 5254.
∗ Electronic address:[email protected] † Electronic address:[email protected]
6
‡ Permanent addres: Centro At´omico Bariloche and Insti- parameters for GaAs: m∗ =0.067me and ǫ=12.4ǫ0.
tuto Balseiro, 8400 S.C. de Bariloche, R´ıo Negro, Ar- 15 A. K. Rajagopal, J. C. Kimball, and M. Banerjee, Phys.
gentina Rev.B 18, 2339 (1978).
1 Fora review, see,e.g., R.M.Dreizler and E.K.U.Gross, 16 E. R¨as¨anen, A. Castro, and E. K. U. Gross Phys. Rev. B
Density functional theory (Springer, Berlin, 1990); U. von 77, 115108 (2008).
Barth, Phys.Scr. T109, 9 (2004). 17 A. D. Becke and M. R. Roussel, Phys. Rev. A 39, 3761
2 Forareview, see,J.P.PerdewandS.Kurth,inA Primer (1989).
in Density Functional Theory (Springer, Berlin, 2003). 18 M. Taut, J. Phys.A 27, 1045 (1994).
3 Forareview,see,e.g., L.P.Kouwenhoven,D.G.Austing, 19 J. B. Krieger, Y. Li, and G. J. Iafrate, Phys. Rev. A 46,
and S. Tarucha, Rep. Prog. Phys. 64 (2001) 701; S. M. 5453 (1992).
Reimann and M. Manninen, Rev. Mod. Phys. 74 (2002) 20 M.A.L.Marques,A.Castro,G.F.BertschandA.Rubio,
1283. Comp.Phys.Comm.151,60(2003);A.Castro,H.Appel,
4 A.K.RajagopalandJ.C.Kimball,Phys.Rev.B15,2819 M. Oliveira, C.A.Rozzi, X.Andrade,F.Lorenzen, M.A.
(1977). L.Marques,E.K.U.Gross,andA.Rubio,Phys.Stat.Sol.
5 B.Tanatar,D.M.Ceperley,Phys.Rev.B39,5005(1989). (b) 243, 2465 (2006).
6 C. Attaccalite, S. Moroni, P. Gori-Giorgi, and G. B. 21 M.Rontani,C.Cavazzoni,D.Bellucci,andG.Goldoni,J.
Bachelet, Phys. Rev.Lett. 88, 256601 (2002). Chem. Phys. 124, 124102 (2006).
7 H. Saarikoski, E. R¨as¨anen, S. Siljam¨aki, A. Harju, M. J. 22 M. Koskinen, M. Manninen, and S. M. Reimann, Phys.
Puska, and R. M. Nieminen, Phys. Rev. B 67, 205327 Rev.Lett. 79, 1389 (1997).
(2003). 23 F. Pederiva, C. J. Umrigar, and E. Lipparini, Phys. Rev.
8 S.Pittalis,E.R¨as¨anen,andM.A.L.Marques,Phys.Rev. B 62, 8120 (2000); ibid68, 089901 (2003).
B 78, 195322 (2008). 24 A.Harju,S.Siljam¨aki,andR.M.Nieminen,Phys.Rev.B
9 S. Pittalis, E. R¨as¨anen, N. Helbig, and E. K. U. Gross, 65, 075309 (2002).
Phys. Rev.B 76, 235314 (2007). 25 N.Helbig,S.Kurth,S.Pittalis, E.R¨as¨anen, andE.K.U.
10 E.R¨as¨anen,S.Pittalis,C.R.Proetto,andE.K.U.Gross, Gross, Phys. Rev.B 77, 245106 (2008).
submitted. 26 This is due to the fact that for N = 1 systems, ψσ(r) =
11 S. Pittalis, E. R¨as¨anen, J. G. Vilhena, and M. A. L. Mar- exp[iθ(r)]pρσ(r),withoutlossofgenerality.Evaluationof
ques, Phys.Rev. A 79, 012503 (2009). the three contributions to Dσ(r) in Eq. (13), for this par-
12 A. D.Becke, J. Chem. Phys. 88, 1053 (1988). ticular case, yields Dσ(r)=0
13 U. von Barth and L. Hedin,J. Phys.C 5, 1629 (1972).
14 We use the effective atomic units (a.u.) with the material
