Kohn-Sham density functional theory for quantum wires in arbitrary correlation regimes Francesc Malet,1 Andr´e Mirtschink,1 Jonas C. Cremon,2 Stephanie M. Reimann,2 and Paola Gori-Giorgi1 1Department of Theoretical Chemistry and Amsterdam Center for Multiscale Modeling, FEW, Vrije Universiteit, De Boelelaan 1083, 1081HV Amsterdam, The Netherlands 2Mathematical Physics, Lund University, LTH, P.O. Box 118, SE-22100 Lund, Sweden (Dated: January 31, 2013) We use the exact strong-interaction limit of the Hohenberg-Kohn energy density functional to construct an approximation for the exchange-correlation term of the Kohn-Sham approach. The resulting exchange-correlation potential is able to capture the features of the strongly-correlated regime without breaking the spin or any other symmetry. In particular, it shows “bumps” (or barriers)thatgiverisetochargelocalizationatlowdensitiesandthatareawell-knownkeyfeature 3 oftheexactKohn-Shampotentialforstrongly-correlatedsystems. Hereweillustratethisapproach 1 for the study of both weakly and strongly correlated model quantum wires, comparing our results 0 2 withthoseobtainedwiththeconfigurationinteractionmethodandwiththeusualKohn-Shamlocal density approximation. n a J I. INTRODUCTION the consequent need of considering larger Hilbert spaces 0 in the calculations. Other wave-function methods like 3 In semiconductor nanostructures, the regime of strong Quantum Monte Carlo7,18,19 (QMC) and density matrix ] correlation is reached when the electronic density be- renormalizationgroup(DMRG),20whichrelytosomeex- el comes low enough so that the Coulomb repulsion be- tent on various approximations, can treat systems larger - comesdominantwithrespecttothekineticenergyofthe than the CI approach, but are still computationally ex- tr electrons. From the purely fundamental point of view, pensive and limited to N (cid:46)100. s the study of the strongly-interacting limit in such sys- ThemuchcheaperKohn-Sham(KS)densityfunctional . at tems is interesting since charge localization, reminiscent theory (DFT),21,22 which allows to treat thousands of m oftheWignercrystallization1 ofthebulkelectrongas, is electrons,isthemethodofchoicetotreatlargerquantum expected to occur at low densities. systems. However,allthecurrentlyavailableapproxima- - d A lot of previous theoretical work on Wigner localiza- tions for the exchange-correlation functional fail to de- n tion in nanostructures has focused on finite-sized quan- scribe the strongly-correlated regime7,23–27 even at the o tumdots(see,forexample,Refs. 2–7),andthecrossover qualitative level. Allowing spin- and spatial-symmetry c fromliquidtolocalizedstatesinthetransportproperties breaking may yield reasonable total energies, without, [ ofthenanostructurehasbeenaddressed.8,9 Inquasione- however, capturing the physics of charge localization in 1 dimensional nanosystems, signatures of Wigner localiza- non-magnetic systems. Moreover, broken symmetry so- v tion were observed experimentally in one-dimensional lutions often yield a wrong characterization of various 3 cleaved-edge overgrowth structures,10 or in the trans- properties and the rigorous KS DFT framework is par- 2 portpropertiesofInSbnanowirequantum-dotsystems.11 tially lost (see, e.g., Refs. 20,23,26). 3 7 More recent experimental work clearly identified the for- KSDFTis,inprinciple,anexacttheorythatshouldbe . mation of Wigner molecules in a one-dimensional quan- abletoyieldtheexactenergyanddensityeveninthecase 1 tum dot that was capacitively coupled to an atomic ofstrongelectroniccorrelation,withoutartificiallybreak- 0 3 force microscope probe.12 Wigner localization has also inganysymmetry. However,whendealingwithpractical 1 been investigated in other 1D systems such as carbon KSDFT,onecouldexpectthatthenon-interactingrefer- : nanotubes.13–15 (For a review, see Ref. 16). Finally, encesystemintroducedbyKohnandShammightnotbe v regarding practical applications, Wigner-localized sys- the best choice when trying to address systems in which i X tems have been shown to be potentially useful,e.g., for the electron-electron interactions play a dominant role. r quantum-computing purposes.13,17 For many years, huge efforts have been made in order to a When trying to model electronic strongly-correlated trytogetabettercharacterizationandunderstandingof systems, however, the commonly employed methodolo- the properties of the exact Kohn-Sham reference system gies encounter serious difficulties of different nature. On (see e.g. Refs. 20,25,28–44). All these works reflected the one hand, the configuration interaction (CI) ap- the large difficulties encountered when trying to obtain proach, despite being in principle capable of describing adequate approximations to describe strong correlation anycorrelationregime,isinpracticelimitedtothestudy in the exact KS theory.45 of small systems with only very few particles due to its An alternative density-functional framework, based highcomputationalcost,whichscalesexponentiallywith on the study of the strongly-interacting limit of the the number of particles, N. Such numerical difficulties Hohenberg-Kohn density functional, was presented in get even worse in the very strongly-correlated limit due Ref.46. Inthisapproach,areferencesystemwithinfinite to the degeneracy of the different quantum states and correlation between the electrons was considered instead 2 of the non-interacting one of Kohn and Sham. The two theperformancesofKSSCEwiththe“exact”CIresults, formalisms can therefore be seen as complementary to with the standard KS LDA method, and discussing KS each other and, indeed, the first results obtained with SCEwithasimplelocalcorrection. Finally,inSec.Vwe this so-called strictly-correlated-electrons (SCE) DFT, draw some conclusions, as well as an outlook for future presently limited to either 1D or spherically-symmetric works. systems, showed its ability to describe systems in the Hartree (effective) atomic units are used throughout extreme strongly-correlated regime with a much better the paper. accuracy than standard KS DFT.46,47 On the downside, however,SCEDFTrequiresthatoneknowsapriorithat thesystemisinthestrong-interactionregime,anditfails II. THEORY AND METHODOLOGY as soon as the fermionic nature of the electrons plays a significantrole.47 Furthermore,theformalismlackssome A. KS and SCE DFT of the appealing properties of the Kohn-Sham approach, such as its capability to predict (at least in principle) In the formulation of Hohenberg and Kohn21 the exact ionization energies. Also, crucial concepts widely ground-state density and energy of a many-electron sys- employed in solid state physics and in chemistry, such as temareobtainedbyminimizingtheenergydensityfunc- the Kohn-Sham orbitals and orbital energies, are totally tional absent in SCE DFT. (cid:90) Very recently, a new approach that combines the E[ρ]=F[ρ]+ drv (r)ρ(r) (1) ext advantages of the KS and the SCE DFT formalisms, consisting in approximating the Kohn-Sham exchange- with respect to the density ρ(r). In Eq. (1) v (r) is correlation energy functional with the strong-interaction ext the external potential and F[ρ] is a universal functional limit of the Hohenberg-Kohn energy density functional, has been proposed.48 Pilot tests of this new “KS SCE” of the density, defined as the minimum of the internal energy (kinetic energy Tˆ plus electron-electron repulsion framework showed that it is able to capture the features of both the weakly and the strongly-correlated regimes Vˆee) with respect to all the fermionic wave functions Ψ in semiconductor quantum wires, as well as the so-called that yield the density ρ(r),49 2k → 4k crossover occurring in between them, while F F keeping (at least for 1D systems) a computational cost F[ρ]= min(cid:104)Ψ|Tˆ+Vˆee|Ψ(cid:105). (2) Ψ→ρ comparable to the one of standard KS DFT with the local-density approximation (LDA). In other words, the In order to capture the fermionic nature of the elec- SCE functional yields a highly non-local approximation tronic density, Kohn and Sham22 introduced the func- for the exchange-correlation energy functional, which is tional T [ρ] by minimizing the expectation value of Tˆ s able to capture key features of strong correlation within alone over all the fermionic wave functions yielding the the KS scheme, without any artificial symmetry break- given ρ(r),49 ing. Themainpurposeofthisworkistofurtherinvestigate T [ρ]= min(cid:104)Ψ|Tˆ|Ψ(cid:105), (3) s this new KS SCE method, by discussing its exact formal Ψ→ρ properties and, for the prototypical case of (quasi)-1D thus introducing a reference system of non-interacting quantum wires, by also performing full CI calculations electronswiththesamedensityasthephysical, interact- to compare electronic densities, total energies and one- ing,one. TheremainingpartofF[ρ],definingtheHartree electron removal energies in different regimes of correla- and the exchange-correlation functionals, F[ρ]−T [ρ]≡ s tion. We find that the KS SCE results are qualitatively E [ρ] ≡ E [ρ] + E [ρ], is then approximated. The Hxc H xc right at all correlation regimes, representing an impor- minimization of the total energy functional E[ρ] with re- tant advance for KS DFT. However, while one-electron spect to the density yields the well-known single-particle removal energies are quite accurate, total energies and Kohn-Sham equations22 ground-state densities are still quantitatively not always satisfactory, and therefore we also discuss the construc- (cid:18) 1 (cid:19) − ∇2+v [ρ](r) φ (r)=ε φ (r), (4) tion of corrections to KS SCE. In particular, we investi- 2 KS i i i gatehereasimplelocalcorrection,which,however,turns out to give rather disappointing results, suggesting that where v (r) ≡ v (r)+δE [ρ]/δρ(r) ≡ v [ρ](r)+ KS ext Hxc ext to further improve KS SCE we need semi-local or fully v [ρ](r)+v [ρ](r)istheone-bodylocalKohn-Shampo- H xc non-local density functionals. tential, with v [ρ](r) and v [ρ](r) being, respectively, H xc The paper is organized as follows. In the next Sec. II the Hartree and the exchange-correlation parts. The so- we describe the KS SCE approach, illustrating and dis- lutions φ of Eqs. (4) are the so-called Kohn-Sham or- i cussingitsfeaturesbeyondwhatwasreportedinRef.48. bitals, which yield the electronic density through the re- InSec.IIIweintroducethequasi-1Dsystemswehavead- lationρ(r)=(cid:80) |φ (r)|2,withthesumrunningonlyover i i dressed,andinSec.IVwepresentourresults,comparing occupied orbitals. Notice that here we work with the 3 original, spin-restricted, KS scheme, in which we have invariance of ρ under the coordinate transformation s→ the same KS potential for spin-up and spin-down elec- f (s), i.e., i trons. The HK functional of Eq. (2) and the KS functional ρ(f (s))df (s)=ρ(s)ds, (9) i i of Eq. (3) can be seen as the particular values at λ = 1 and at λ=0 of a more general functional F [ρ] in which or, equivalently, that the probability of finding the elec- λ the coupling-strength interaction is rescaled with a real tron i at fi(s) is equal to that of finding the electron parameter λ, i.e., “1” at s. At the same time, the fi(s) must satisfy group properties that ensure the indistinguishability of the N F [ρ]= min(cid:104)Ψ|Tˆ+λVˆ |Ψ(cid:105). (5) electrons.55,57 λ ee Ψ→ρ The functional VSCE[ρ] can then be written in terms ee of the co-motion functions f as55,59 A well-known exact formula for the Hartree-exchange- i correlation functional E [ρ] is50,51 Hxc (cid:90) ρ(s)N(cid:88)−1 (cid:88)N 1 VSCE[ρ] = ds (cid:90) 1 (cid:90) 1 ee N |f (s)−f (s)| EHxc[ρ]= (cid:104)Ψλ[ρ]|Vee|Ψλ[ρ](cid:105)dλ≡ Veλe[ρ]dλ, (6) i=1 j=i+1 i j 0 0 1(cid:90) (cid:88)N 1 = dsρ(s) , (10) where Ψλ[ρ] is the minimizing wave function in Eq. (5). 2 |s−fi(s)| In the strictly-correlated-electrons DFT (SCE DFT) i=2 formalism, one considers the strong-interaction limit of just as T [ρ] is written in terms of the Kohn-Sham or- s the Hohenberg-Kohn functional, λ → ∞, which corre- bitals φ (r). The equivalence of the two expressions for sponds to the functional52–55 VSCE[ρ]iin Eq. (10) has been proven in Ref. 59. ee Since in the SCE system the position of one electron VSCE[ρ]≡ min(cid:104)Ψ|Vˆ |Ψ(cid:105), (7) ee ee determines all the other N − 1 relative positions, the Ψ→ρ net repulsion felt by an electron at position r due to i.e., theminimumoftheelectronicinteractionaloneover the other N − 1 electrons becomes a function of r it- all the wave functions yielding the given density ρ(r). self. For a given density ρ (r), this effect can be ex- 0 This limit has been first studied in the seminal work actly transformed46,55,57 into a local one-body effective of Seidl and coworkers52–54, and later formalized and external potential v [ρ ](r) that compensates the to- SCE 0 evaluated exactly in a rigorous mathematical way in talCoulombforceoneachelectronwhenalltheparticles Refs.47,55–57. ThefunctionalVSCE[ρ]alsodefinesaref- areattheirrespectivepositionsf [ρ ](r),i.e.,suchthat55 ee i 0 erencesystemcomplementarytothenon-interactingone oinffitnhieteKlyo-hcnor-Srehlaamtedkienleetcitcroennes,rgwyiTths[ρz]e,roonkeincoemticpoesneedrgbyy. ∇v [ρ ](r)=(cid:88)N r−fi[ρ0](r) . (11) SCE 0 |r−f [ρ ](r)|3 Thisimpliesthat,analogouslyasinasetofconfinedclas- i=2 i 0 sicalrepulsivecharges,whicharrangethemselvesseeking for the stable spatial configuration that minimizes their Intermsoftheclassical-chargeanalogue,vSCE[ρ0](r)can interaction energy, in the SCE reference system the po- thus be seen as an external potential for which the total sitionofoneelectronuniquelydeterminesthepositionof classical potential energy theremainingones,alwaysundertheconstraintimposed N−1 N N byEq. (7)thatthedensityateachpointisequaltothat (cid:88) (cid:88) 1 (cid:88) E (r ,...r )≡ + v [ρ ](r ) of the quantum-mechanical system with λ=1, ρ(r). pot 1 N |r −r | SCE 0 i i j Moreprecisely,thefunctionalVSCE[ρ]isconstructed55 i=1 j=i+1 i=1 ee (12) by considering that the admissible configurations is minimum when the electronic positions reside on the of N electrons in d dimensions are restricted to subset R , i.e., when r = f [ρ ](r) or, equivalently, a d−dimensional subspace Ω0 of the full classical when theΩa0ssociated densiity ati e0ach point is equal to Nd−dimensional configuration space. A generic point ρ (r). Foranarbitrarydensityρ(r),thepotential-energy of Ω has the form 0 0 density functional defined as RΩ0(s)=(f1(s),.....,fN(s)), (8) (cid:90) ESCE[ρ]≡VSCE[ρ]+ v [ρ ](r)ρ(r)dr (13) pot ee SCE 0 where s is a d-dimensional vector that determines the position of, say, electron “1”, and f (s) (i = 1,...,N), i will then satisfy the stationarity condition with f1(s) = s, are the so-called co-motion functions, δESCE[ρ]/δρ(r)(cid:12)(cid:12) =0, i.e., we will have that which determine the position of the i-th electron as a pot ρ=ρ0 function of s. The co-motion functions are implicit non- sloocluatliofunnoctfiaonsaeltsooffdtihffeergeinvteinaldeeqnusaittyioρn(srt)h,4a6t,5e5,n5s7u,5r8eatnhde δVδeSρeC(rE)[ρ](cid:12)(cid:12)(cid:12)(cid:12) =−vSCE[ρ0](r). (14) ρ=ρ0 4 Notice that Eq. (14) involves the functional derivative of therefore allowing one to address both the weakly- and a highly non-local implicit functional of the density, de- the strongly-interacting regime, as well as the crossover finedbyEqs.(9)-(10). This,howeverturnsouttoreduce between them.48 Indeed, from the scaling properties60 of toalocalone-bodypotentialthatcanbeeasilycalculated thefunctionalsF[ρ],T [ρ]andVSCE[ρ]itderivesthatthe s ee fromtheintegrationofEq.(11)oncetheco-motionfunc- approximation of Eq. (18) becomes accurate both in the tionsareobtainedviaEq.(9). Thisshortcuttocompute weak-andinthestrong-interactionlimits,whileprobably thefunctionalderivativeofVSCE[ρ]isextremelypowerful lesspreciseinbetween. Tousethescalingrelations60 one ee for including strong-correlation in the KS formalism.48 defines, for electrons in D dimensions, a scaled density ρ (r)≡γDρ(γr) γ >0. γ B. Zeroth-order KS-SCE approach We then have47,60 T [ρ ] = γ2T [ρ] (19) Equations(11)and(14)showhowtheeffectsofstrong s γ s correlation, captured by the limit λ → ∞ of Fλ[ρ] and VeSeCE[ργ] = γVeSeCE[ρ] (20) rigorously represented by the highly non-local functional F[ρ ] = γ2F [ρ], (21) γ 1/γ VSCE[ρ],areexactlytransferredintotheone-bodypoten- ee tial vSCE[ρ]. The KS SCE approach to zeroth order con- where F1/γ[ρ] means60 that the Coulomb coupling con- sists in using this property to approximate the Hartree- stantλinFλ[ρ]ofEq.(5)hasbeensetequalto1/γ. We exchange-correlation term of the Kohn-Sham potential then see that both sides of Eq. (18) tend to Ts[ργ] when as γ → ∞ (high-density or weak-interaction limit) and to VSCE[ρ ] when γ →0 (low-density or strong-interaction ee γ δEHxc[ρ] ≈v˜ [ρ](r), v˜ [ρ](r)≡−v [ρ](r). limit). δρ(r) SCE SCE SCE Standard KS DFT emphasizes the non-interacting (15) shellstructure,properlydescribedthroughthefunctional Noticethatwehavedefinedv˜ [ρ](r)=−v [ρ](r),as T [ρ], but it misses the features of strong correlation. SCE SCE s here we seek an effective potential for KS theory, which SCE DFT, on the contrary, is biased towards localized corresponds to the net electron-electron repulsion acting “Wigner-like” structures in the density, accurately de- on an electron at position r, while the effective poten- scribed by VSCE[ρ], missing the fermionic shell struc- ee tial for the SCE system of Eq. (11) compensates the net ture. Manyinterestingsystemslieinbetweentheweakly electron-electron repulsion. and the strongly interacting limits, and their complex More rigorously, by considering the λ→∞ expansion behavior arises precisely from the competition between of the integrand of Eq. (6) one obtains52–55,58 the fermionic structure embodied in the kinetic energy andcorrelationeffectsduetotheelectron-electronrepul- VZPE[ρ] sion. By implementing the exact v˜ [ρ](r) potential in Vλ→∞[ρ]=VSCE[ρ]+ ee√ +O(λ−p), (16) SCE ee ee λ the Kohn-Sham scheme, we thus let these two factors compete in a self-consistent procedure.48 where the acronym “ZPE” stands for “zero-point en- One should also notice that while the KS SCE ap- ergy”, and p ≥5/4 – see Ref. 58 for further details. By proachdoesnotuseexplicitlytheHartreefunctional,the inserting the expansion of Eq. (16) into Eq. (6) one ob- correct electrostatics is still captured, since VSCE[ρ] is ee tains an approximation for EHxc[ρ], the classical electrostatic minimum in the given density ρ. Moreover, the potential v˜ [ρ](r) stems from a wave SCE EHxc[ρ]≈VeSeCE[ρ]+2VeZePE[ρ]+... (17) function (the SCE one55,58) and is therefore completely self-interaction free. We consider here only the first term, corresponding to a Finally, another neat property of the zeroth-order KS zeroth-order expansion around λ = ∞, i.e., E [ρ] ≈ Hxc SCE approach is that it always yields a lower bound to VSCE[ρ], which yields Eq. (15) for the corresponding ee theexactground-stateenergyE0 =E[ρ0],whereρ0isthe functional derivatives. exact ground-state density. In fact, for any given ρ the Taking into account the definition of the functional right-hand side of Eq. (18) is always less or equal than VSCE[ρ], Eq. (7), the zeroth-order KS SCE is equivalent ee the left-hand side, as the minimum of a sum is always toapproximatetheminimizationoverΨintheHKfunc- larger than the sum of the minima. As a consequence, tional of Eq. (2) as for ρ=ρ we have the inequality 0 (cid:90) (cid:90) min(cid:104)Ψ|Tˆ+Vˆee|Ψ(cid:105) ≈ min(cid:104)Ψ|Tˆ|Ψ(cid:105)+ min(cid:104)Ψ|Vˆee|Ψ(cid:105) E[ρ ]=F[ρ ]+ ρ v ≥T [ρ ]+VSCE[ρ ]+ ρ v , Ψ→ρ Ψ→ρ Ψ→ρ 0 0 0 ext s 0 ee 0 0 ext = T [ρ]+VSCE[ρ]. (18) (22) s ee which becomes even stronger when ones minimizes the TheKSSCEapproachthustreatsboththekineticenergy functionalontheright-hand-sidewithrespecttotheden- and the electron-electron repulsion on the same footing, sity within the self-consistent zeroth-order KS SCE pro- combiningtheadvantagesofbothKSandSCEDFTand cedure. It should be noted that this property implies an 5 importantdifferencewithrespecttothevariationalwave- III. MODEL AND DETAILS OF THE function methods (such as HF, CI, QMC and DMRG), CALCULATIONS which, instead, provide an upper bound to the exact ground-state energy. We consider N electrons in the quasi-one-dimensional (Q1D) model quantum wire of Refs. 27,62, C. Local correction to zeroth-order KS SCE Hˆ =−1(cid:88)N ∂2 +N(cid:88)−1 (cid:88)N w (|x −x |)+(cid:88)N v (x ), 2 ∂x2 b i j ext i As preliminary found in Ref. 48 and further shown i=1 i i=1 j=i+1 i=1 (30) in Sec. IV, the zeroth-order KS SCE yields results that in which the effective electron-electron interaction is ob- are qualitatively correct in the strong-correlation regime tained by integrating the Coulomb repulsion on the lat- (representing a significative conceptual advance for KS eral degrees of freedom,62 and is given by DFT), but still with quantitative errors, which become √ smaller and smaller as correlation increases. An impor- π (cid:18) x2 (cid:19) (cid:16) x (cid:17) tant issue is thus to add corrections to Eq. (18). One wb(x)= 2b exp 4b2 erfc 2b . (31) can, more generally, decompose F[ρ] as F[ρ]=T [ρ]+VSCE[ρ]+T [ρ]+Vd[ρ], (23) The parameter b fixes the thickness of the wire, set to s ee c ee b=0.1 throughout this study, and erfc(x) is the comple- where Tc[ρ] (kinetic correlation energy) is mentaryerrorfunction. Theinteractionwb(x)hasalong- T [ρ]=(cid:104)Ψ[ρ]|Tˆ|Ψ[ρ](cid:105)−T [ρ], (24) range coulombic tail, wb(x→∞)=1/x, and is finite at c s the origin, where it has a cusp. As in Ref. 27, we con- i.e., the difference between the true kinetic energy and sider an external harmonic confinement v (x)= 1ω2x2 ext 2 the Kohn-Sham one, and Vd[ρ] (electron-electron decor- in the direction of motion of the electrons. The wire ee relation energy) is can be characterized by an effective confinement-length parameter L such that Vd[ρ]=(cid:104)Ψ[ρ]|Vˆ |Ψ[ρ](cid:105)−VSCE[ρ], (25) ee ee ee 4 1 i.e., the difference between the true expectation of Vˆee ω = L2, vext(x)= 2ω2x2. and the SCE value. A “first-order” approximation for T [ρ]+Vd[ρ] can be obtained from Eq. (17), c ee T [ρ]+Vd[ρ]≈2VZPE[ρ], (26) A. Zeroth-order KS SCE c ee ee and can be, in principle, included exactly using the for- The co-motion functions f (x) can be constructed by i malism developed in Ref. 58, but other approximations, integratingEqs.(9)foragivendensityρ(x),52,56,57choos- e.g. in the spirit of Ref. 61, can also be constructed. ing boundary conditions that make the density between Here we consider an even simpler approximation, two adjacent strictly-correlated positions always inte- Tc[ρ] + Vede[ρ] ≈ ELC[ρ], where ELC[ρ] is a local term grate to 1 (total suppression of fluctuations),52 thatincludes,ateachpointofspacer,thecorresponding correction for a uniform electron gas with the same local (cid:90) fi+1(x) density ρ(r), i.e., ρ(x(cid:48))dx(cid:48) =1, (32) (cid:90) fi(x) ELC[ρ]= ρ(r)(cid:2)tc(ρ(r))+vede(ρ(r))(cid:3)dr. (27) and ensuring that the fi(x) satisfy the required group properties.52,55,57 This yields In Eq. (27) t (ρ) and vd (ρ) are the kinetic correlation c ee (cid:40) energyandtheelectron-electrondecorrelationenergyper N−1[N (x)+i−1] x≤a f (x)= e e N+1−i particle of an electron gas with uniform density ρ, corre- i N−1[N (x)+i−1−N] x>a , e e N+1−i sponding to (33) tc(ρ)+vede(ρ)=(cid:15)xc(ρ)−(cid:15)SCE(ρ), (28) where the function Ne(x) is defined as where (cid:15) (ρ) is the usual electron-gas exchange- (cid:90) x xc N (x)= ρ(x(cid:48))dx(cid:48), (34) correlation energy and (cid:15) (ρ) is the indirect part (ex- e SCE pectation of Vˆ minus the Hartree energy) of the SCE −∞ ee interaction energy per electron of the uniform electron and ak =Ne−1(k). Equation (11) becomes in this case gas with density ρ. This correction makes the approxi- N mate internal energy functional (cid:88) v˜(cid:48) [ρ](x)= w(cid:48)(|x−f (x)|)sgn(x−f (x)). (35) SCE b i i F[ρ]=Ts[ρ]+VeSeCE[ρ]+ELC[ρ] (29) i=2 becomeexactinthelimitofuniformdensity,similarlyto Wethensolveself-consistentlytheKohn-Shamequations what the LDA functional does in standard KS DFT. (4) with the KS potential v (x)=v (x)+v˜ [ρ](x), KS ext SCE 6 wherev˜ [ρ](x)isobtainedbyintegratingEq.(35)with SCE 0 the boundary condition v˜ [ρ](|x| → ∞) = 0. As said, SCE -0.5 we work in the spin-restricted KS framework, in which (cid:161) each spatial orbital is doubly occupied. gy -1 er -1.5 n E e -2 (cid:161)xc B. The configuration interaction method (CI) C -2.5 fit (cid:161) S SCE d -3 data (cid:161)SCE n In the configuration interaction calculations, the full c a -3.5 many-body wavefunction is expanded as a linear combi- x -4 nation of Slater determinants, constructed with the non- -4.5 interacting harmonic oscillator orbitals. A matrix repre- 0 2 4 6 8 10 12 14 sentation of the Hamiltonian in this basis is then numer- Wigner-Seitz radius r ically diagonalized to find the eigenstates of the system. s The number of possible ways to place N particles in a given set of orbitals increases rapidly as a function of FIG. 1: (color online) The indirect SCE energy (cid:15) (r ) for N, such that only small particle numbers are tractable. SCE s the 1D gas [interaction of Eq. (31) and b=0.1] is compared Also, the stronger the interaction, the more basis or- totheparametrizedQMCdata66fortheexchange-correlation bitals are generally required to obtain a good approxi- energy(cid:15) (r ). FortheSCEenergyweshowbothournumer- xc s mation. For the present physical system, about 20–40 ical results and the fitting function of Eqs. (38)-(39). orbitals were needed to get converged solutions, which resulted in Hilbert space dimensions in the range 105– 106. For a more detailed description of the method, see have then evaluated the limit N →∞ at fixed density ρ e.g. Refs. 63,64. to obtain the bulk value. The details of this calculation are reported in Appendix A. In Fig. 1 we show our numerical results for b = 0.1 C. KS LDA compared to the parametrized66 QMC results for the exchange-correlationenergy(cid:15) (r )ofEqs.(36)-(37). We xc s We have performed Kohn-Sham LDA calculations us- see that, as it should be, (cid:15)SCE(rs)≤(cid:15)xc(rs) everywhere. ing the exchange-correlation energy per particle (cid:15)xc = For large rs we find that the SCE data are very close (cid:15) + (cid:15) for a 1D homogeneous electron gas with the to the QMC parametrization, with differences of the or- x c renormalized Coulomb interaction wb(x), as detailed in der of ∼ 0.1%. Notice also√that at rs = 0 we have Ref. 27. The exchange term (cid:15)x is given by (cid:15)SCE(0)=(cid:15)xc(0)=(cid:15)x(0)=− 4bπ. This is due to the fact that in the r → 0 limit at fixed b the first-order per- s 1(cid:90) +∞ dq turbation to the non-interacting gas is just a constant, (cid:15) (r )= v (q)[S (q)−1] , (36) x s 2 2π b 0 so that every normalized wave-function yields the same −∞ result for the leading term. We have parametrized our wherevb(q)istheFouriertransformoftheinteractionpo- data for (cid:15)SCE(rs) as tential, S (q) is the non-interacting static structure fac- 0 tor, and r ≡ 1 .65 To increase the numerical stability, (cid:15)SCE(ρ)=ρq(2bρ), (38) s 2ρ we have interpolated between the Taylor expansions of with (cid:15) (r ) at small and large r up to order 14. For the cor- rexlatsion term we have usedsthe results of Casula et al.,66 q(x)=A ln(cid:18) a1x+a2x2 (cid:19), (39) who have parametrized their QMC data as 1 1+a x+a x2 3 2 (cid:15) (r )=− rs ln(1+Dr +Erγ2), (37) and A1 = 0.9924534, a√2 = 1.55176743, a3 = c s A+Brsγ1 +Crs2 s s 2.025166778, a1 = a3 − a22A1π. This fit is valid for all values of b, since the scaling of Eq. (38) is exact for the where the different parameters are given in Table IV of SCE energy. The fitting function is also shown in Fig. 1 Ref. 66 for several values of b. for the case b=0.1. D. KS SCE with local correction IV. RESULTS WehaveobtainedtheindirectSCEenergyperelectron Figure 2 shows the electron densities for N = 4 and (cid:15)SCE(ρ)neededinEq.(28)byfirstcomputingtheindirect different effective confinement lengths L = 2ω−1/2 ob- (cid:15)drop(ρ,N) for a 1D droplet with N electrons, uniform tained with the KS SCE, the CI and the KS LDA ap- SCE densityρandradiusR= N,asdescribedinRef.56. We proaches. One can see that the three methods show 2ρ 7 qualitative agreement in the weakly-correlated regime, represented here in panel (a) by the case L = 1. The densities have N/2 peaks, given by the Friedel-like oscil- lations with wave number 2keff, where keff =πρ˜/2 is the 1.2 F F effective Fermi wavenumber, determined by the average density in the bulk of the trap ρ˜. 0.8 L=1 KKSS -LSDCEA As the confinement length of the wire increases, the CI interactions start to become dominant and, whereas the 0.4 N=4 KS SCE and the CI results are still in qualitative agree- ment,theLDAclearlyprovidesaphysicalwrongdescrip- 0 tion of the system. Indeed, one can see from panel (b) that whereas the densities obtained from the KS SCE 0.6 L=15 and the CI methods develop a four-peak structure, cor- 2 respondingtochargelocalizationandindicatingthatthe L/0.4 system enters the crossover between the weakly and the ρ strongly correlated regimes (the 2k → 4k crossover), 0.2 F F the KS LDA yields a flat density. This is a typical er- 0 ror of local and semilocal density functionals that shows L=70 up also in bond breaking (yielding wrong molecular dis- 0.6 sociation curves) and in systems close to the Mott in- sulating regime. In such cases, better total energies are 0.4 obtained by using spin-dependent functionals and allow- ing symmetry breaking. This, however, does not yield a 0.2 satisfactoryphysicaldescriptionofsuchsystems,missing many key features and giving a wrong characterization 0-8 -4 0 4 8 of several properties (see, e.g., Refs. 20,23,26). 2x/L When the system becomes even more strongly- correlated, hererepresentedbyL=70, theKSSCEgets closertotheCIresult,withdensitiesthatclearlypresent FIG. 2: (color online) Electron densities for N = 4 and N peaks, corresponding to charge localization. The KS L = 1, 15 and 70, obtained with the KS SCE, CI and LDA LDA density is now very delocalized and almost flat in approaches. The results are given in units of the effective the scale of Fig. 2. In order to obtain charge localization confinement length L=2ω−1/2. within the restricted KS scheme, the self-consistent KS potential must build “bumps” (or barriers) between the electrons. These barriers are a very non-local effect and are known to be a key property of the exact Kohn-Sham In particular, Vieira44 has shown that the exact potential, as discussed in Refs. 38 and 29. exchange-correlation potential for a 1D Hubbard chain In Fig. 3 we show that the self-consistent KS SCE with hopping parameter t and on-site interaction U, ob- scheme builds, indeed, the above-mentioned barriers in tained by inversion from the exact many-body solution, the corresponding Kohn-Sham potentials, which we plot always oscillates with frequency 4k , while the density F together with the corresponding densities for N =4 and oscillations undergo a 2k → 4k crossover with in- F F N = 5 for L = 70. One can see that each of the N creasing U/t. The crossover in the density is thus due peaks in the density corresponds to a minimum in the to the increase in the amplitude of the oscillations of the KS potential, which is separated from the neighboring xc potential. In Fig. 4 we show the KS SCE exchange- ones by barriers or “bumps”, at whose maxima the KS correlation potentials for N = 4 electrons in the weakly potential has a discontinuous (but finite57) first deriva- (L=2)andstrongly(L=70)correlatedregimes. Wesee tive. ThenumberofsuchbarriersisthusequaltoN−1, thattheKSSCEself-consistentresultsareinqualitative andtheybecomemorepronouncedwithincreasingcorre- agreementwiththefindingsofVieira:44theoscillationsin lation, enhancing the corresponding charge localization. the xc potential have essentially a frequency 4k also in F Notice that the discontinuous first derivative of the KS theweakly-correlatedcase,withamplitudethatincreases SCE potential at the barrier maxima is a feature due to withincreasingL[duetothescalingofEqs.(19)-(21)the the classical nature of the SCE potential, and it is not parameterLplaysherearolesimilartoU/tfortheHub- expected to appear in the exact KS potential (indeed, it bard chain]. In the two lower panels of the same figure does not appear in any of the available calculations of we also further clarify the 2k → 4k crossover in the F F the “exact” KS potential obtained by inversion). KSframework: weseethatthe4k regimeinthedensity F It is also interesting to make a connection between oscillations occurs when the barriers in the total KS po- our results and the recent work on the KS exchange- tential (due to the large oscillations of the xc potential) correlation potential for the 1D Hubbard chains.44,67,68 are large enough to create classically-forbidden regions 8 0 2 4 6 82 N=4 L=2 L=70 N=4 0.4 E 80 SCxc v L=70 0.2 78 2 2 2/ L/ρ 0 104 LKS L=2 N=5 v L=70 0.4 S ρ K L=70 v 0.2 100 0 0 2 4 6 2x/L x FIG. 3: (color online) Self-consistent Kohn-Sham potentials obtained with the KS-SCE method for N = 4 and N = 5, with effective confinement length L = 70 (blue solid lines). FIG. 4: (color online) Top panel: the self-consistent KS SCE Thecorrespondingdensitiesarealsoshown(reddottedlines). exchange-correlation (xc) potential for N =4 at weak corre- Notice that for the sake of clarity only the results for x > 0 lation (L=2) and strong correlation (L=70). In the inset, are shown. The results are given in units of the effective the oscillating part of the xc potential at L = 2 is zoomed confinement length L=2ω−1/2. in. Middle panel: the total self-consistent KS SCE potential (blue,solidline),thecorrespondingdensity(reddottedline), andthetwooccupiedKSeigenvalues(greendashedhorizontal lines) for the weakly-correlated L = 2 wire. In this case, we inside the trap for the occupied KS orbitals. see that in the KS system there are no classically-forbidden In Table I we report the total energies obtained with regions inside the trap. Bottom panel: the same as in the the three approaches, KS SCE, CI and KS LDA, for dif- middle panel for the strongly-correlated L=70 wire. In this ferent values of the parameters L and N. It can be seen case,weclearlyseetheclassically-forbiddenregionsinsidethe thatintheweakly-correlatedregime,representedhereby trap created by the barriers in the KS SCE potential. The results are given in arbitrary units. L=1 and 2, the error made by the KS SCE approach is larger than the one corresponding to the KS LDA. The results also clearly show that, as previously discussed, KS SCE is always a lower bound to the total energy. As N L KS SCE CI KS LDA thesystembecomesmorecorrelated,theresultsobtained 2 2 1.81 2.49 2.59 with the KS SCE and the CI approaches become closer 2 15 0.0942 0.106 0.130 toeachother, whereasthevaluegivenbytheKSLDAis 2 70 0.0112 0.0115 0.0182 less accurate, as one could have inferred from the corre- 4 1 25.08 28.42 28.57 sponding densities shown in panels b) and c) of Fig. 2. 4 2 8.46 10.60 10.68 4 15 0.491 0.541 0.580 In the exact Kohn-Sham theory, the highest occu- 4 70 0.0602 0.0629 0.0771 pied KS eigenvalue is equal to minus the exact chem- 5 15 0.787 0.871 0.915 ical potential from the electron-deficient side,69,70 i.e., 5 70 0.099 0.102 0.121 µ− =E −E . InTableIIwecomparethehighestoc- N−1 N cupiedKSeigenvalueobtainedwiththeKSSCEandthe TABLE I: Comparison of the total energies obtained with KS LDA approaches with the values of EN −EN−1 cal- theKSSCE,CIandKSLDAapproachesfordifferentvalues culated from the total energies given by the CI method, of the particle number N and effective-confinement length corresponding to the same values of N and L given in L=2ω−1/2. Table I. One can see that in this case the KS SCE gives 9 N L KS SCE CI KS LDA 2 2 1.65 1.99 2.56 2 15 0.104 0.097 0.263 2 70 0.0126 0.0111 0.04087 a) 4 1 11.26 11.86 12.56 4 2 4.08 4.65 5.02 N=8 4 15 0.248 0.256 0.453 L=70 4 70 0.0318 0.0304 0.06909 5 15 0.325 0.330 0.539 5 70 0.0408 0.0391 0.08172 b) TABLE II: For the same systems of Table I, we compare the N=8 highest occupied KS eigenvalues obtained from KS SCE and L=150 KS LDA with the full CI values of E −E . N N−1 S ρ K c) v N=16 good results also in the weakly-correlated regime. In L=150 the strongly-correlated limit, the KS SCE and the CI re- sults show an agreement similar to that observed in the corresponding total energies. KS LDA, as usual, yields d) too high eigenvalues, due to the too fast decay of the N=32 exchange-correlation potential for |x|→∞. L=150 As mentioned earlier, the numerical cost of the CI method increases exponentially with the number of par- ticles, and this limitation becomes stronger as the cor- x relations become dominant. In the calculations reported above, for the 5-electron case with L =70 we diagonal- ized a matrix where the eigenvectors had a dimension of about 3.5×105. While it is technically possible to treat FIG.5: (coloronline)ElectrondensityandcorrespondingKS larger matrices, the rapid growth of the basis size still SCE potential for different particle numbers N and effective efficiently limits the number of particles one can han- confinementlengthsL. AsinFig.3,onlytheresultsforx>0 dle. (For N =6 electrons, using the same basis orbitals, are shown. The results are given in arbitrary units. the corresponding dimension is roughly 2.6×106.) The KS SCE method, on the contrary, has a numerical cost (in 1D) comparable to the one of KS LDA, therefore al- lowing to study strongly-correlated systems with much larger particle numbers. In Fig. 5 we show the electron densities and corresponding KS potentials obtained with the KS SCE method for N =8, 16, and 32, for different 0.6 KS-SCE KS-SCE+LDA values of L: in panels a) and b) we see how, at fixed CI KS LDA number of particles N =8, the bumps in the KS poten- N=2 L=20 tial and the amplitude of the density oscillations become 0.4 larger with increasing L. For fixed effective confinement 2 length L=150, we see from panels b), c) and d) how in- L/(cid:108) creasing the particle number N leads to less pronounced features of strong correlation, according to the scaling of 0.2 Eqs. (19)-(21). Finally, we have tested the local correction to the zeroth-orderKSSCEdiscussedinSecs.IICandIIID:as 0 -6 -4 -2 0 2 4 6 we see in the case N =2 and L=20 reported in Fig. 6, 2x/L theresultsfortheself-consistentdensitiesareverydisap- pointing,layinginbetweentheKSSCEandthestandard KS LDA values. This is due to the fact that, similarly FIG. 6: (color online) Electron density for the case N = 2 and L=20. The “exact” CI result is compared with the KS to the standard KS LDA case, this simple local correc- LDA, the KS SCE and the KS SCE with local correction of tion cannot capture the physics of the intermediate and Secs. IIC and IIID (KS SCE+LDA) results. The results are strong-correlation regime, so that its inclusion worsens giveninunitsoftheeffectiveconfinementlengthL=2ω−1/2. the results of KS SCE. In future work we will explore semi-local and fully non-local corrections to KS SCE. 10 V. CONCLUSIONS AND PERSPECTIVES tion for Scientific Research (NWO) through a Vidi grant andbytheSwedishResearchCouncilandtheNanometer We have used the exact strong-interaction limit of Structure Consortium at Lund University (nmC@LU). the Hohenberg-Kohn functional to approximate the exchange-correlationenergyandpotentialofKohn-Sham DFT. By means of this so-called KS SCE approach, we Appendix A: SCE for the uniform Q1D electron gas have addressed quasi-one-dimensional quantum wires in theweak,intermediateandstrongregimeofcorrelations, FollowingRef.56,wehavecomputedtheSCEindirect comparing the results with those obtained by using the Coulomb interaction energy per electron (cid:15)drop(ρ,N) of configuration interaction and the KS local density ap- SCE a 1D droplet of uniform density ρ and radius R = N, proximation. In the weakly-correlated regime, the three 2ρ where N is the number of electrons, approaches give qualitatively similar results, with elec- tronic densities showing N/2 peaks, associated with the (cid:18) (cid:19) 2 ρ 2bρ double occupancy of the single-particle levels that domi- (cid:15)drop(ρ,N)= ρv˜SCE(2ρb,N)− u , (A1) SCE N ee π 1 N nate the system. In this regime, KS LDA performs over- all better than KS SCE. As correlations become domi- where nant, the KS SCE and the CI densities start to develop additional maxima, corresponding to charge-density lo- N (cid:18) (cid:19) calization, whereas the KS LDA provides a qualitatively v˜SCE(x,N)= π (cid:88)(N −i)ei2/x2erfc i (A2) wrongdescriptionofthesystem,yieldingaveryflat,delo- ee 2x x i=1 calized, density. Wehavealsoinvestigatedasimplelocal correction to KS SCE, which, however, gives very dis- is the rescaled SCE energy of the droplet56 and the sec- appointing results. In future works we will thus explore ondtermintheright-hand-sideofEq.(A1)isitsHartree semi-local and fully non-local corrections to KS SCE. energy, with The Kohn-Sham potential of the KS SCE approach shows “bumps” that are responsible for the charge local- (cid:90) ∞(cid:18)sink(cid:19)2 ization and are a well-known feature of the exact Kohn- u1(x)= k ek2x2E1(k2x2)dk, (A3) Sham potential of strongly-correlated systems. More- 0 over,theassociatedKSSCEexchange-correlationpoten- and tialshowstherightasymptoticbehaviour,sinceitisself- interaction free as it is constructed from a wave function (cid:90) ∞ e−tx (the SCE one48,55,58). This way, KS SCE is able to also E1(x)= t dt. give rather accurate chemical potentials. Notice that, as 1 shownbystudiesofone-dimensionalHubbardchains,the Sincethefunctionu (x)isnumericallyunstable,wehave 1 2k →4k crossover in the density is a very challenging F F interpolated between its small-x expansion through or- taskforKSDFTfornon-magneticsystems.44,67 Thefact ders O(x5), that KS SCE is able to capture this crossover is thus a very remarkable and promising feature. πx4 πx2 1 1 (cid:18)3(cid:19) Crucial for future applications is calculating VSCE[ρ] u<(x)= − + π3/2x−πlog(x)− πψ(0) , ee 1 16 4 2 2 2 and v˜ [ρ](r) also for general two- and three- SCE dimensionalsystems. Anenticingroutetowardsthisgoal with ψ(0)(cid:0)3(cid:1) ≈ 0.036489974, and its large-x expansion involvesthemass-transportation-theoryreformulationof 2 the SCE functional,57 in which VSCE[ρ] is given by the through orders O(x−16), ee maximum of the Kantorovich dual problem, π3/2 64π π3/2 (cid:90) (cid:88)N N(cid:88)−1(cid:88)N 1  u>1(x)= 1209600x15 − 14189175x14 + 131040x13 − max u(r)ρ(r)dr : u(r )≤ , u  i |ri−rj| 16π π3/2 16π π3/2 2π i=1 i=1 j>i + − + − 405405x12 15840x11 51975x10 2160x9 945x8 whereu(r)=v˜SCE[ρ](r)+C,andC isaconstant.57 This π3/2 4π π3/2 π π3/2 π π3/2 is a maximization under linear constraints that yields in + − + − + − + , 336x7 315x6 60x5 15x4 12x3 3x2 2x one shot the functional and its functional derivative. Al- though the number of linear constraints is infinite, this switching between them at x=0.584756. formulation may lead to approximate but accurate ap- We have then evaluated numerically the limit N →∞ proaches to the construction of VSCE[ρ] and v˜ [ρ](r), ee SCE of Eq. (A1) at fixed ρ. We have found that the conver- as very recently shown by Mendl and Lin.71 gence is reasonably fast: for example, taking N = 105 Acknowledgments yields results with a relative accuracy of 10−6. Our nu- mericalresultshavebeenfittedwiththefunctionq(x)of WethankM.Seidlforinspiringdiscussions. Thiswork Eqs. (38)-(39). was financially supported by the Netherlands Organiza-