Density functional description of Coulomb blockade: Adiabatic or dynamic exchange-correlation?
DDensity functional description of Coulomb blockade:Adiabatic or dynamic exchange-correlation?
Zhen-Fei Liu and Kieron Burke Molecular Foundry and Materials Sciences Division,Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA Departments of Chemistry and Physics, University of California, Irvine, California 92697, USA (Dated: October 16, 2018)Above the Kondo temperature, the Kohn-Sham zero-bias conductance of an Anderson junctionhas been shown to completely miss the Coulomb blockade. Within a standard model for the spectralfunction, we deduce a parameterization for both the onsite exchange-correlation potential and thebias drop as a function of the site occupation that applies for all correlation strengths. We use ourresults to sow doubt on the common interpretation of such corrections as arising from dynamicalexchange-correlation contributions.
I. INTRODUCTION
Electron transport [1] through molecular junctions isusually formulated using non-equilibrium Green’s func-tion (NEGF)[2] or an equivalent scattering formalism [3].In atomistic simulations of electron transport [4, 5], den-sity functional theory (DFT) [6, 7] is often used in com-bination with NEGF, producing the standard model[8]in which the interacting one-body Green’s function andcoupling matrices are replaced by their Kohn-Sham (KS)counterparts [4, 5]. Whether or not ground-state KS-DFT can in principle, under certain conditions, producethe exact conductance remains open[8], even at zero tem-perature and in the linear response regime. Thus thereare two sources of error which can be difficult to dis-entangle: (i) replacing the interacting one-body Green’sfunction with the KS Green’s function, and (ii) replacingthe exact KS Green’s function by one from an approxi-mate functional[9].Model Hamiltonians provide excellent testbeds for thispurpose, especially when highly accurate results areknown and can serve as benchmarks for larger KS-DFTcalculations. The Anderson impurity model [10] has at-tracted much attention recently[11–18]. Due to its gener-ality, broad applicability, and exact solvability, studyingtransport through an Anderson junction is useful in un-derstanding errors in atomistic calculations and in con-structing accurate functionals for electron transport.The KS system is defined as a junction with U = 0 (seebelow), but whose onsite energy is chosen to make its siteoccupation match that of the interacting junction, whichcan be found exactly using the Bethe ansatz[19]. Atzero temperature and in the linear response regime, theKS conductance, i.e., the conductance of non-interactingelectrons in the KS potential, is exact[11, 12, 14], i.e., thefirst error is zero. This can be attributed to the Friedel-Langreth sum rule [20, 21], which implies the transmis-sion is a simple function of occupation. By reverse en-gineering of the exact solution[12], one can study whichfeatures must be present in any approximation in orderto generate an accurate transmission[13]. For example,the derivative discontinuity[22] is crucial as correlations grow, and how rapidly it is approached as U grows isdetermined by the charge susceptibility of the system atparticle-hole symmetry [13]. This derivative discontinu-ity has also been shown to be related to more generalblockade phenomena [23, 24].However, all this changes above the Kondotemperature[25], when the Kondo transmission peak isnegligible, and the sum rule no longer applies. Accuratesolutions are available [1], and the spectral functionexhibits two Coulomb peaks in the strongly correlatedregime. Between the two peaks, the linear-responseconductance is very low due to Coulomb blockade.This temperature regime better mimics what happensin real molecular junctions than the zero-temperaturelimit where a Kondo plateau dominates. But the low-temperature KS conductance always satisfies the sumrule even when the physical conductance does not, andso is qualitatively incorrect between the two Coulombpeaks. Having spent years figuring out how to get theKondo plateau into DFT, the difficulty is now to get ridof it at finite temperature [26].It has long been argued[27–29] that time-dependentDFT produces dynamic corrections to the KS conduc-tance, and that these are necessary to produce an accu-rate transmission. Making this identification, Ref. [15]shows that a simple model for the exchange-correlation(XC) kernel of time-dependent density functional theory(TDDFT) reproduces the right corrections, at least inthe strongly correlated limit. This was further shown todepend, somewhat surprisingly, on the onsite occupationalone.Here, we argue that the corrections that turn theKS conductance into the true conductance of an An-derson junction can plausibly be considered as a non-local static correction that in principle could be extractedfrom ground-state DFT. This possibility was suggestedlong ago [28] as the only alternative to dynamic TDDFTeffects for altering a resonance in the transmission ofsuch a junction. The origin of such an effect, withinground-state DFT, is the well-documented counteractingXC field that is significant for certain molecules and solidsin response to a long-range electric field. The origin of a r X i v : . [ c ond - m a t . s t r- e l ] M a y this field is the localization of orbitals on specific sitesand the appearance of step-like features in the inducedexchange-correlation potential that (correctly) reduce thepolarization, relative to that of standard local, semilocal,and hybrid functionals[30]. These steps are reasonablyaccurately given by Hartree-Fock or optimized effectivepotential (OEP) calculations, because of their explicitorbital dependence[31]. They are in fact a field-inducedderivative discontinuity.To make this argument, we first accurately parametrizethe onsite XC potential for temperatures that are low,but at which there is no Kondo plateau in the con-ductance, i.e., in the Coulomb blockade regime. Ourparametrization is designed to work for all correlationstrengths, not just for strong correlation. Next, we ac-curately parametrize the XC bias drop that ensures thata static KS calculation reproduces the physical conduc-tance, again including both weak and strong correlationregimes. Finally, we explain how this can be interpretedin terms of the known counteracting XC fields. II. ANDERSON JUNCTION AT LOWTEMPERATURE LIMIT
To begin, the Hamiltonian for the Andersonjunction[10] consists of an interacting impurity site cou-pled to identical featureless left and right leads: H = ε ˆ n + U ˆ n ↑ ˆ n ↓ + H leads + H T (1)with ε the on-site energy and U the Coulomb repulsionwhen the site is doubly occupied, ˆ n = ˆ n ↑ + ˆ n ↓ . Here H leads is the Hamiltonian of the two leads, and H T is theircoupling to the site. In the wide band limit, the effectof tunneling is incorporated into an energy-independentconstant Γ [32]. There are two dimensionless parame-ters: ( µ − ε ) /U , with µ the chemical potential of theleads, which moves the system on- and off-resonance, and u = U/ Γ, which switches the system between weakly andstrongly correlated. Below the Kondo temperature T K ,a characteristic temperature dependent on u , the exactspectral function has two Coulomb peaks and one Kondopeak whose width is related to T K [33], and the sum ruleapplies. Above T K , the Kondo peak disappears, the sumrule is violated, and the conductance comes from the twoCoulomb peaks only [33].Although no exact solution is available above T K , theGreen’s function on the central impurity site can be ac-curately approximated as [1] G ( ω ) = (cid:88) i =1 , n i ω − ε i + i Γ / , (2)where ε = ε, ε = ε + U, n = 1 − n/ , and n = n/ n ≡ (cid:104) ˆ n C (cid:105) the occupation on the impurity site, andthe coupling to the leads causes the broadening of Γ / ansatz forthe exact solution, as it captures the right physics and accurately mimics the numerical exact solution above T K [1]. We also introduce a specific low temperature limit,namely β − →
0, but β − (cid:29) k B T K , where β is the in-verse temperature and k B the Boltzmann constant. Be-cause T K depends exponentially on the parameters, it istypically much smaller than any other temperature scale.We are in a regime above T K , but at temperatures farsmaller than any of the other energy scales of the prob-lem. For example, Kurth and Stefanucci [15] discusseda case where U = 10 , Γ = 1 , and temperature τ = 0 . T K = 0 . τ = 0 . T K , to facilitate comparison of our ap-proach (static correction) and that in Ref. [15] (dynamiccorrection). Ref. [15] studied these effects for a rangeof temperatures but only for strong correlation; here weare interested in all correlation strengths, and find thatlow temperatures are well approximated by a specific low-temperature limit. In this limit, we find accurate analyticparametrizations with relative ease.The spectral function A ( ω ) = − G ( ω ) depends onthe occupation n on the impurity site. Therefore n mustbe determined self-consistently: n = 1 π (cid:90) µ −∞ dωA ( ω ) = π − θ π + θ − θ , (3)where tan θ i = 2( ε i − µ ) / Γ and, in the low temperaturelimit, levels are fully occupied up to µ . Thanks to thislimit, Eq. (3) is a closed-form for n in terms of the pa-rameters. The linear-response transmission is then [15]: T = − Γ2 (cid:90) ∞−∞ dω π dfdω A ( ω ) β − → −−−−−→ Γ4 π A ( ω = µ ) . (4)In Fig. 1, n and T are plotted as a function of µ , forseveral values of u . For U (cid:29) Γ, a plateau in n develops,i.e., a Coulomb blockade, which corresponds to the lowconductance region between the two Coulomb peaks at µ ≈ (cid:15) and µ ≈ (cid:15) + U . This is very different from theKondo regime, where there is always a Kondo plateau inconductance, as discussed in Ref. [12]. III. KOHN-SHAM ANDERSON JUNCTION
We now construct a KS Anderson junction, i.e., onewith U = 0, that generates the same occupation n as theinteracting system by replacing ε with ε S [ n ] = ε + v HXC [ n ],and analyze its conductance. Because the occupation onthe impurity site is a number, the functional dependenceis simply a function of n . The KS Green’s function hasonly one pole: G S ( ω ) = 1 ω − ε S [ n ] + i Γ / . (5) T n ( µ - ε )/Uu=0.1u=1u=10 FIG. 1. Upper panel: Conductance T as a function of ( µ − ε ) /U , plotted in unit of G = 2 e /h = 1 /π , the conductancequantum. lower panel: n as a function of ( µ − ε ) /U . Different u = U/ Γ are used.
The Hartree-exchange-correlation (HXC) potential v HXC [ n ] is defined such that: n = 1 π (cid:90) µ −∞ dωA S ( ω ) , (6)where A S ( ω ) = − G S ( ω ) is the spectral function ofthe KS system. By definition[34], the KS occupationmatches the physical one, i.e., the left hand sides of Eqs.(3) and (6) are identical, for a given set of µ, ε, U, and Γ.Applying this condition yields: v HXC = µ − ε − Γ2 tan (cid:104) π n − (cid:105) . (7)Note that this is insufficient to define v HXC [ n ] because ofthe presence of µ − ε [If Eq. (3) could be inverted to find µ − ε as an explicit function of n , as is trivial numerically,this would suffice]. In Ref. [15], reverse engineering tofind v HXC [ n ] is done at three different temperatures above T K , but for a fixed large u (=10). Here the low tem-perature limit simplifies the algebra, and we study the u -dependence explicitly. Following the ideas of Ref. [13],we parametrize v HXC [ n ] as v app HXC [ n ] = U (cid:18) π tan − [ σ ( u )( n − (cid:19) , (8)where σ ( u ) is a parameter, determined by the ex-act condition of charge susceptibility at particle-holesymmetry[35]: χ = U ∂n∂µ (cid:12)(cid:12)(cid:12)(cid:12) n =1 = 4 u (1 + u )( π + 2 tan − u ) , (9)where Eq. (3) was used. The derivative is taken atthe particle-hole symmetry point, n = 1, or equivalently where (cid:15) = µ − U/
2. Imposing this condition fixes theparameter σ in Eq. (8) as σ ( u ) = π/χ − π / (4 u ). Thisprocedure could be simply generalized to finite temper-ature, by using the finite-temperature susceptibility, butthis is not the main concern of the present work. Webriefly show and discuss finite-temperature effects in Sec-tion VI using the numerical solution, without analyticallyparametrizing the temperature-dependent functional.
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 v HX C [ n ] / U n u=1 u=2u=5 u=20 FIG. 2. v HXC [ n ] as a function of n for various u = U/ Γ. Solidlines are results from Eq. (7) with n calculated from Eq. (3),dashed lines are the approximation of Eq. (8). In Fig. 2, both the v HXC [ n ] that precisely reproducesthe onsite occupation determined from Eq. (3) and theparametrization of Eq. (8) are plotted as a function of n , for different u . For U (cid:29) Γ, a step in v HXC devel-ops, reflecting the onset of the derivative discontinuityas u → ∞ [34]. These kinds of features were discussedin Refs. [12, 13], at zero temperature. In the low tem-perature limit discussed in this work, these rules still ap-ply. An important point is that, unlike Ref. [15], ourparametrization applies for all values of U (both weakand strong correlation), and appears to fail only whenthe site is almost entirely empty or doubly-occupied. Ourresults match those of Ref. [15] when the temperature islow and the correlation is strong.The KS conductance has the same form as in Eq. (4),but with A replaced by A S . We compare the two quan-tities in Fig. 3. Here the KS conductance always satis-fies the Friedel-Langreth sum rule, i.e., T S = sin ( πn/ U (cid:29) Γ due to the ansatzof Eq. (2). For U (cid:28) Γ, the KS conductance is quite accu-rate, but for larger U , it is wildly inaccurate. Whenevercorrelation is significant (or, equivalently in this model,the coupling to the leads is weak), the KS conductance isa large overestimate of the physical conductance. In par- T ( µ - ε )/U u=0.1u=1u=10 FIG. 3. Physical conductance (solid lines) from Eq. (4) andKS conductance (dashed lines) for various values of u = U/ Γ,plotted in unit of G = 2 e /h = 1 /π , the conductance quan-tum. ticular, it misses entirely the Coulomb blockade effect,producing no drop at all between the Coulomb peaks.All this was beautifully shown in Ref. [15] for strong cor-relation for several temperatures. Although the Friedel-Langreth sum rule does not apply at finite temperatures,the discrepancy remains qualitatively the same. IV. STATIC CORRECTION TO KOHN-SHAMCONDUCTANCE
It is no surprise that the KS conductance is not cor-rect under such conditions, as it is extracted entirely fromequilibrium DFT[36]. One can formally identify the cor-rections to the KS conductance[8], and capture them interms of an XC correction to the applied bias: δI = T δV = T S δV S , (10)where δV is an applied bias, δV S = δV + δV XC is the biasexperienced by the KS system, and δI is the induced cur-rent in both. The formal result of Ref. [12] can be statedthat δV XC = 0 due to the sum rule at zero temperaturewhen the Kondo peak is present, and the pioneering workof Ref. [15] interpreted the XC corrections in δV S as dy-namical corrections arising from TDDFT, i.e., correctionsthat would not appear in any static DFT calculation.Here we suggest that a static correction from ground-state DFT can explain the results of Kurth andStefanucci[15] which they interpreted solely in terms ofTDDFT. We propose the following approximation for anon-local XC bias: δV app S δV = n [1 − cos(2 ζ )] + 2 n / (cid:2) u + cot ζ ) (cid:3) ( n π ) (11) for n <
1, where ζ = n ( π + 2 n tan − u/n ) and n and n are swapped if n >
1. This specific choice will bejustified shortly. For the present, we simply plot it inFig. 4. For weak correlation at any filling, or for n or n < .
5, this KS bias is close to the applied bias. Butfor n near 1 and U (cid:29) Γ, δV S is noticeably smaller than δV . δ V S / δ V n u=0.1u=0.5u=1u=5u=20 FIG. 4. δV app S as a function of n for various correlationstrengths u = U/ Γ. Solid lines are “exact”, as defined byEq. (12). Dashed lines are approximate, as in Eq. (11).
Using δV app S in Eq. (11) with n calculated self-consistently from Eq. (8), we now plot the approximateconductance in Fig. 5. We see that it gives us essen-tially the correct conductance of the Anderson junction,including the severe corrections needed to reproduce theCoulomb blockade. -2 -1 0 1 2 3 T ( µ - ε )/U u=0.1u=1 u=10 FIG. 5. Solid lines: physical conductances from Eq. (4).Dashed lines: conductances using approximations for boththe local HXC potential Eq. (8) and non-local XC bias Eq.(11).
To understand how this was achieved, i.e., to derive δV app S , we first ask what XC drop would exactly repro-duce the conductance implied by Eq. (4). From Eq. (10),we see δV S A S ( ω = µ ) = δV A ( ω = µ ) . (12)The δV S /δV is an implicit functional of n via ( µ − ε ) /U .For example, at particle-hole symmetry, we find explicitly δV S δV (cid:12)(cid:12)(cid:12)(cid:12) n =1 = 11 + u , β − → . (13)This “exact” δV S is plotted in Fig. 4 as a function of n .The KS bias is close to the applied bias when the siteis nearly empty or doubly-occupied. But near n = 1, as U increases, δV S /δV drops rapidly from 1 to 0, i.e., theapplied bias is almost entirely screened by the Coulombblockade.To use δV XC to correct the KS conductance, one needsto express it as an explicit functional of the density n .The conductance is proportional to A ( ω = µ ), whichdepends explicitly on ε , the external onsite potential (thestrength of electron-electron interaction, u = U/ Γ, entersinto the formula as a parameter). We start with Eq. (3),and propose the following approximation: θ ( n ) ≈ ( n −
12 ) π + 2 n tan − u, (14)for n <
1, and n and n are swapped if otherwise. ThenEq. (3) can be inverted to find ε [ n ] = µ + Γ cot ζ/ ε [ n ] into Eqs. (2) and (4), we findour expression for δV S /δV . One can verify that Eq. (11)satisfies the relation Eq. (13). The approximation usedin Ref. [15] corresponds to θ ≈ π/ u limit, whereasEq. (14) works wherever δV XC is significant, even forsmall u or weak correlation [In the language of Ref. [15], R = 2 n tan − u/ ( πn )]. V. ORIGIN OF NON-LOCAL CORRECTION
To understand how the XC bias drop might be ex-tracted from a ground-state DFT calculation, we reviewthe origin of Eq. (10). Start with Kubo response of thesystem to the external field, and here we follow the no-tation in Refs. [28] and [37]: δ j ( r ; ω ) = (cid:90) d r (cid:48) ˆ σ irr ( r , r (cid:48) ; ω ) δ E tot ( r (cid:48) ; ω ) (15)where the left-hand side is the current density in responseto a perturbing external electric field at frequency ω , andˆ σ irr is the irreducible, frequency-dependent, nonlocal con-ductivity tensor of the many-electron problem. Since,in time-dependent current DFT (TDCDFT), the time-dependent KS system must produce the same current- density response, we also have δ j ( r ; ω ) = (cid:90) d r (cid:48) ˆ σ S ( r , r (cid:48) ; ω ) [ δ E tot ( r (cid:48) ; ω ) + δ E XC ( r (cid:48) ; ω )] , (16)where δ E tot contains both the external and Hartree elec-tric fields. Here ˆ σ S is the KS conductivity tensor and δ E XC is the exchange-correlation contribution to the fieldfelt by the KS system.The current is defined as an integral of the current den-sity over a cross section: δI ( z ; ω ) = (cid:82) S dS δ j ( r ; ω ), where z is the direction along the current flow. Great care withthe order of limits must be taken when applying this for-mula to our problem[37]. The response formula is true forany ω , but we wish to deduce the steady-state current.Thus ω is kept finite, but reduced to 0 at the end of thecalculation. In that limit, the nonlocal conductivity andthe current become coordinate-independent[28, 37, 38],and are just the transmission T in Eq. (10). Integra-tion of the fields over a large region including the de-vice just yields the net voltage drop across the device, δV = lim L →∞ (cid:82) L − L dz (cid:82) S dS δ E ( r ; ω ). Thus we recoverEq. (10).Now comes the tricky part. For any local (or semilo-cal) approximation to XC, if the density deep inside theleads is symmetric far from the molecule, the XC fieldsmust be also symmetric, so that there can be no net XCbias drop. This is the case in all standard DFT calcula-tions of transport[8, 28, 39]. But Ref. [28] argues thatthere are two possible sources of δV XC : either highly non-local ground-state effects ( ω = 0) or dynamical TDDFTcorrections which fail to vanish as ω → a priori reason to dis-count a non-local static contribution, which has provedmuch more straightforward (if expensive) to implementfor both the optical response of insulators and polariz-abilities of long-chain polymers. Thus it is important tounderstand if the effects seen in Ref. [15] can be under-stood in this manner. If so, then they can be searchedfor in orbital-dependent calculations of transport, suchas GW, etc.We note here that the crucial element for producingthis counteracting XC field is not whether or not the sys-tem is a metal or insulator (traditional definitions basedon bulk conductance[51] break down for nanoscale sys-tems), but rather whether or not charge is localized ona site. This is exactly what is being modeled in the An-derson junction. Furthermore, while TDCDFT is neededto derive Eq. (10), the establishment of the steady cur-rent is done by the time-dependence of the KS equationsthemselves, without any need for dynamic XC contribu-tions. As shown by Ref. [37], a steady current is es-tablished within even a simple Hartree calculation, i.e.,one including no XC effects at all. This shows that thebasic feature of a steady current occurs with the time-dependent KS equations in the absence of any dynamicalXC contributions. -1-0.8-0.6-0.4-0.2 0 0.2 -1 -0.5 0 0.5 1 1.5 2 δ V X C / δ V ( µ - ε )/Uu=0.1u=0.5u=1u=5u=20 FIG. 6. δV XC as a function of µ for various correlationstrengths u = U/ Γ. Solid lines are “exact” as defined byEq. (12) and dashed lines are approximations using bothEqs. (11) and (8).
To illustrate how this interpretation makes sense, wealso plot δV XC in Fig. 6 as a function of µ − ε . Forlarge U the XC bias almost completely cancels the ap-plied field, just as happens for well-separated moleculesin response to applied fields[30]. Because we have per-formed our analysis for all U , we see that this remainstrue at all correlation strengths, for all ε < µ < ε + U , i.e.,in the region between the Coulomb peaks. Thus there isalways an opposing field in this region, which can be un- derstood as the origin of Coulomb blockade in terms ofdensity functionals. Finally, we note that outside this re-gion, we see regions where the XC bias has the same signas the applied bias. We suspect this is an artifact of thesimplicity of the model. Either the Anderson junctionitself is oversimplified so its results cannot be trusted inthis region, or possibly this is the non-applicability of oursimple model for the spectral function [Eq. (2)] in thisregime. In either case, this effect should be treated withcaution unless it is also seen in a real-space calculation. VI. FINITE TEMPERATURE EFFECTS
The primary focus of our work has been the low-temperature limit, because of its relevance to standardelectronic structure calculations. However, it is relativelystraightforward to repeat our calculations at higher tem-peratures, as was done in Ref. [15]. To do that, Eq. (3)needs to be modified as: n = 2 (cid:90) ∞−∞ dω π f ( ω ) A ( ω ) , (17)where f ( ω ) is the Fermi function: f ( ω ) = 1 / (1+exp[( ω − µ ) /τ ]), with τ being the temperature. Also, the finitetemperature version of Eq. (4) needs to be used, i.e., T = − Γ2 (cid:90) ∞−∞ dω π dfdω A ( ω ) . (18)At finite temperature, we have no simple closed-formparametrization for the density or HXC potential andeverything we show in this section is numerical. For easeof comparison, we use the same temperatures as those ofRef. [15], namely τ = 0 . , . , and 1.0, respectively, with u = 10. We also show our low-temperature limit result.In Fig. 7, we show conductance and occupation ofthe impurity site as a function of ( µ − ε ) /U for differenttemperatures (similar results can be found in Fig. 1 ofRef. [15]). In the lower panel, we see that for tempera-tures below about 0.2, our low-temperature limit resultis indistinguishable from the numerical result. On theother hand, the conductance shows substantially greatersensitivity, both exactly and at the KS level, due to thepresence of the Fermi function in Eq. (18). The Friedel-Langreth sum rule applies only in the low-temperaturelimit. However, the figure also shows that our low-temperature result is approached as the temperature isreduced.Finally, in Fig. 8 we show the XC drop δV XC /δV asa function of n , for different temperatures. One can seethat large temperatures wash out the sharp effects of thederivative discontinuity in the strongly correlated limit,as expected. The XC drop is less sensitive to temperaturethan the conductance. Overall, the qualitative featuresof the XC counterfield remain, but are weakened by in-creasing temperature, as the derivative discontinuity isrounded off. T n ( µ - ε )/Ulow- ττ =0.1 τ =0.2 τ =1.0 FIG. 7. Upper panel: conductance T as a function of ( µ − ε ) /U , plotted in unit of G , for different temperatures τ at afixed u = U/ Γ = 10. Solid lines are “exact” and dashed linesare (uncorrected) KS conductance. The black line is for thelow-temperature limit and is the same as the red line in Fig. 1.Lower panel: n as a function of ( µ − ε ) /U . By definition, theKS occupation always matches that of the exact occupationat all temperature. -1-0.8-0.6-0.4-0.2 0 0.2 0 0.5 1 1.5 2 δ V X C / δ V n low- ττ =0.1 τ =0.2 τ =1.0 FIG. 8. Finite temperature generalization of Fig. 4 (but with δV XC plotted rather than δV S ): δV XC as a function of n , fordifferent temperatures τ at a fixed u = U/ Γ = 10. Black lineis the low-temperature limit result.
VII. CONCLUDING REMARKS
We conclude with a discussion of the relevance of theseresults for DFT calculations of transport. First we notethat we have provided accurate parametrizations for boththe onsite potential and the XC bias for the junctionat low (but above T K ) temperatures, that apply for alljunctions, not just those strongly coupled to the leads. These are available for any future calculations of impu-rity models of transport. If the parameters (Γ, U , (cid:15) ) canbe extracted from standard ground-state DFT calcula-tions, the results could be compared with the accompa-nying DFT transport calculations to identify limitationsof the pure DFT approximations. If U/ Γ is large, thestandard DFT approach will often yield a large overesti-mate of the conductance, as is often seen in comparisonwith experiment[8].But, possibly more importantly, we have suggestedthat an alternative origin of the corrections to the KSconductance is non-local static XC effects. This is veryimportant for understanding the limitations of standardDFT calculations of conductance and finding ways toimprove them. Our reading suggests that, rather thanlooking to TDCDFT for corrections to the conductance,one could look instead at orbital-dependent ground-state approximations[30, 47, 52–54]. Such calculations canyield two important improvements in calculations in thisarea: (i) the correct positioning of orbital energies in themolecular region relative to the leads, as is already well-known[55], and (ii) a finite XC bias drop that correctsthe KS conductance.Unfortunately, due to the structureless nature of theleads in the Anderson model, there is no way at presentto distinguish between the interpretation of the XC cor-rections presented in Ref. [15] and those presented here.In a more realistic model, the charge is depleted fromthe end of one lead and is increased at the end of theother, producing a net dipole. Such a dipole contributesto the net bias drop in the Hartree potential. But, withan orbital-dependent functional, a counteracting effectoccurs in the exchange bias. This can be seen mostlyclearly when the molecule is only weakly coupled to theleads. Then the form of this counteracting field is stepsthat are proportional to the applied field but that inhibitcharge from moving from the molecule to the leads. Suchsteps are well-known in the absence of an applied field,between two distinct species, where they ensure dissocia-tion into charge neutral fragments, and are consequencesof the derivative discontinuity [22]. Our steps are field-induced versions of the same thing. Only a more detailedcalculation, such as those of Refs. [56] and [17], can yieldsuch information. But at least we have provided a reasonto look for such corrections when non-local ground-stateapproximations are being applied.To conclude, we have shown that, in the case of the An-derson junction, simple parametrizations of both the on-site HXC potential and the XC bias in linear response canreproduce the low-temperature Coulomb blockade, forall correlation strengths, not just the strongly correlatedregime discussed in Ref. [15]. Furthermore, the XC biascan be interpreted as a non-local XC response to an ap-plied electric field, and is not necessarily a dynamic TD-CDFT effect. Qualitatively similar effects should occurin more realistic descriptions of molecular transport[55].Our approximate formulas could prove useful in othercontexts, or as a check on XC approximations applied totransport problems.
VIII. ACKNOWLEDGEMENT
We thank Justin Smith, Klaus Capelle, and FerdinandEvers for helpful discussions. Work at Berkeley Lab wassupported by the U.S. Department of Energy, Office of Basic Energy Sciences, Materials Sciences and Engineer-ing Division, under Contract [1] H. J. W. Haug and A.-P. Jauho,
Quantum Kineticsin Transport and Optics of Semiconductors , 2nd ed.(Springer, Berlin, 2007).[2] Y. Meir and N. S. Wingreen, Phys. Rev. Lett. , 2512(1992).[3] D. S. Fisher and P. A. Lee, Phys. Rev. B , 6851 (1981).[4] J. Taylor, H. Guo, and J. Wang, Phys. Rev. B , 245407(2001).[5] M. Brandbyge, J.-L. Mozos, P. Ordej´on, J. Taylor, andK. Stokbro, Phys. Rev. B , 165401 (2002).[6] W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965).[7] K. Burke, J. Chem. Phys. , 150901 (2012).[8] M. Koentopp, C. Chang, K. Burke, and R. Car, J. Phys.:Condens. Matter , 083203 (2008).[9] C. Toher, A. Filippetti, S. Sanvito, and K. Burke, Phys.Rev. Lett. , 146402 (2005).[10] P. W. Anderson, Phys. Rev. , 41 (1961).[11] G. Stefanucci and S. Kurth, Phys. Rev. Lett. , 216401(2011).[12] J. P. Bergfield, Z.-F. Liu, K. Burke, and C. A. Stafford,Phys. Rev. Lett. , 066801 (2012).[13] Z.-F. Liu, J. P. Bergfield, K. Burke, and C. A. Stafford,Phys. Rev. B , 155117 (2012).[14] P. Tr¨oster, P. Schmitteckert, and F. Evers, Phys. Rev.B , 115409 (2012).[15] S. Kurth and G. Stefanucci, Phys. Rev. Lett. , 030601(2013).[16] G. Stefanucci and S. Kurth, Phys. Status Solidi B ,2378 (2013).[17] F. Evers and P. Schmitteckert, Phys. Status Solidi B ,2330 (2013).[18] F. Evers and P. Schmitteckert, Eur. Phys. Lett. ,47012 (2013).[19] P. B. Wiegmann and A. M. Tsvelick, J. Phys. C: SolidState Phys. , 2281 (1983).[20] J. Friedel, Nuovo Cimento Suppl , 287 (1958).[21] D. C. Langreth, Phys. Rev. , 516 (1966).[22] J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz,Phys. Rev. Lett. , 1691 (1982).[23] K. Capelle, M. Borgh, K. K¨arkk¨ainen, and S. M.Reimann, Phys. Rev. Lett. , 010402 (2007).[24] S. Kurth, G. Stefanucci, E. Khosravi, C. Verdozzi, andE. K. U. Gross, Phys. Rev. Lett. , 236801 (2010).[25] M. Pustilnik and L. Glazman, J. Phys.: Condens. Matter , R513 (2004).[26] S. Kurth and G. Stefanucci, private communication. ().[27] N. Sai, M. Zwolak, G. Vignale, and M. Di Ventra, Phys.Rev. Lett. , 186810 (2005).[28] M. Koentopp, K. Burke, and F. Evers, Phys. Rev. B ,121403 (2006). [29] G. Vignale and M. Di Ventra, Phys. Rev. B , 014201(2009).[30] S. K¨ummel, L. Kronik, and J. P. Perdew, Phys. Rev.Lett. , 213002 (2004).[31] S. K¨ummel and L. Kronik, Rev. Mod. Phys. , 3 (2008).[32] In Refs. [12] and [13], we used Γ = Γ L = Γ R . Here, weuse Γ = Γ L + Γ R , consistent with Ref. [15], which is morerelevant to the discussions in this work. ().[33] R. N. Silver, J. E. Gubernatis, D. S. Sivia, and M. Jarrell,Phys. Rev. Lett. , 496 (1990).[34] D. Carrascal, J. Ferrer, J. C. Smith, and K. Burke, “TheHubbard dimer: A density functional case study of amany-body problem,” ArXiv: 1502.02194 (2015).[35] A. C. Hewson, The Kondo problem to heavy fermions (Cambridge University Press, Cambridge, 1997).[36] N. D. Mermin, Phys. Rev. , A1441 (1965).[37] A. Kamenev and W. Kohn, Phys. Rev. B , 155304(2001).[38] H. U. Baranger and A. D. Stone, Phys. Rev. B , 8169(1989).[39] G. Stefanucci, C.-O. Almbladh, S. Kurth, E. K. U. Gross,A. Rubio, R. van Leeuwen, N. E. Dahlen, and U. vonBarth, in Time-Dependent Density Functional Theory ,edited by M. A. L. Marques, C. A. Ullrich, F. Nogueira,A. Rubio, K. Burke, and E. K. U. Gross (Springer, 2006).[40] X. Gonze, Ph. Ghosez, and R. W. Godby, Phys. Rev.Lett. , 4035 (1995).[41] X. Gonze, Ph. Ghosez, and R. W. Godby, Phys. Rev.Lett. , 294 (1997).[42] R. Resta, Rev. Mod. Phys. , 899 (1994).[43] G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. , 601 (2002).[44] B. Champagne, E. A. Perp`ete, S. J. A. van Gisbergen,E.-J. Baerends, J. G. Snijders, C. Soubra-Ghaoui, K. A.Robins, and B. Kirtman, J. Chem. Phys. , 10489(1998).[45] S. J. A. van Gisbergen, P. R. T. Schipper, O. V. Grit-senko, E. J. Baerends, J. G. Snijders, B. Champagne,and B. Kirtman, Phys. Rev. Lett. , 694 (1999).[46] M. van Faassen, P. L. de Boeij, R. van Leeuwen, J. A.Berger, and J. G. Snijders, Phys. Rev. Lett. , 186401(2002).[47] C. D. Pemmaraju, S. Sanvito, and K. Burke, Phys. Rev.B , 121204 (2008).[48] J. A. Berger, P. L. de Boeij, and R. van Leeuwen, Phys.Rev. B , 035116 (2007).[49] J. Jung, P. Bokes, and R. W. Godby, Phys. Rev. Lett. , 259701 (2007).[50] N. T. Maitra, I. Souza, and K. Burke, Phys. Rev. B ,045109 (2003).[51] W. Kohn, Phys. Rev. , A171 (1964). [52] T. K¨orzd¨orfer, M. Mundt, and S. K¨ummel, Phys. Rev.Lett. , 133004 (2008).[53] A. Ruzsinszky, J. P. Perdew, G. I. Csonka, G. E. Scuseria,and O. A. Vydrov, Phys. Rev. A , 060502(R) (2008).[54] A. Ruzsinszky, J. P. Perdew, and G. I. Csonka, Phys.Rev. A , 022513 (2008). [55] S. Y. Quek, L. Venkataraman, H. J. Choi, S. G. Louie,M. S. Hybertsen, and J. B. Neaton, Nano Lett. , 3477(2007).[56] P. Schmitteckert and F. Evers, Phys. Rev. Lett.100