Transport through correlated systems with density functional theory

TL;DR

This paper reviews recent advances in density functional theory (DFT) for quantum transport, especially through strongly correlated systems, and introduces a steady-state DFT framework that captures both Kondo and Coulomb blockade regimes.

Contribution

It develops a steady-state DFT approach with xc potential and bias corrections, unifying descriptions of Kondo and Coulomb blockade phenomena in quantum transport.

Findings

Correctly captures the Kondo plateau at zero temperature.

Highlights the importance of xc bias correction in Coulomb blockade.

Provides an accurate parametrization for the Anderson model.

Abstract

We present recent advances in Density Functional Theory (DFT) for applications to the field of quantum transport, with particular emphasis on transport through strongly correlated systems. We review the foundations of the popular Landauer-B\"uttiker(LB)+DFT approach. This formalism, when using approximations to the exchange-correlation (xc) potential with steps at integer occupation, correctly captures the Kondo plateau in the zero bias conductance at zero temperature but completely fails to capture the transition to the Coulomb blockade (CB) regime as temperature increases. To overcome the limitations of LB+DFT the quantum transport problem is treated from a time-dependent (TD) perspective using TDDFT, an exact framework to deal with nonequilibrium situations. The steady-state limit of TDDFT shows that in addition to an xc potential in the junction, there also exists an xc correction to the applied bias. Open shell molecules in the CB regime provide the most striking examples of the importance of the xc bias correction. Using the Anderson model as guidance we estimate these corrections in the limit of zero bias. For the general case we put forward a steady-state DFT which is based on the one-to-one correspondence between the pair of basic variables steady density on and steady current across the junction and the pair local potential on and bias across the junction. Like TDDFT, this framework also leads to both an xc potential in the junction and an xc correction to the bias. Unlike in TDDFT, these potentials are independent of history. We highlight the universal features of both xc potential and xc bias corrections for junctions in the CB regime and provide an accurate parametrization for the Anderson model at arbitrary temperatures and interaction strengths thus providing a unified DFT description for both Kondo and CB regimes and the transition between them.

Equations

$$\left[-\frac{{\nabla }^2}{2}+v_0(\mathbf{r})+v_{\mathrm{H}}[n](\mathbf{r})+v_{\mathrm{xc}}[n](\mathbf{r})\right]{\psi }_k(\mathbf{r})={\epsilon }_k{\psi }_k(\mathbf{r})$$

$$n(\mathbf{r})=\sum _{k}f({\epsilon }_k)|{\psi }_k(\mathbf{r}){|}^2$$

$$v_{\mathrm{H}}[n](\mathbf{r})=\int {\mathrm{d}}^3{r}^{\prime }w(\mathbf{r}-{\mathbf{r}}^{\prime })n({\mathbf{r}}^{\prime }),$$

$$v_{\mathrm{xc}}[n](\mathbf{r})=\frac{\delta E_{\mathrm{xc}}[n]}{\delta n(\mathbf{r})}.$$

$$\mbox{\boldmath$\mathcal{H}$}_{s}=\left[\begin{array}[]{ccc}{\mathbf{H}}_{LL}&{\mathbf{H}}_{LC}&0\\ {\mathbf{H}}_{CL}&{\mathbf{H}}_{CC}&{\mathbf{H}}_{CR}\\ 0&{\mathbf{H}}_{RC}&{\mathbf{H}}_{RR}\end{array}\right]$$

$${\psi }_k(\mathbf{r})=\sum _{jq}{c}_{k,jq}⟨\mathbf{r}|iq⟩$$

$$n({\bf r})=\sum_{jq}\sum_{j^{\prime}q^{\prime}}\mbox{\boldmath$\rho$}^{\rm eq}_{jq,j^{\prime}q^{\prime}}\langle jq|{\bf r}\rangle\langle{\bf r}|j^{\prime}q^{\prime}\rangle.$$

$$\mbox{\boldmath$\rho$}^{\rm eq}_{jq,j^{\prime}q^{\prime}}=2\sum_{k}f(\varepsilon_{k})c_{k,jq}^{*}c_{k,j^{\prime}q^{\prime}}\;,$$

$$\left((\omega+i\eta)\mathbbm{1}-\mbox{\boldmath$\mathcal{H}$}_{s}\right)\mbox{\boldmath$\mathcal{G}$}(\omega)=\mathbbm{1}$$

$$\mbox{\boldmath$\mathcal{G}$}(\omega)=\left[\begin{array}[]{ccc}{\mathbf{G}}_{LL}(\omega)&{\mathbf{G}}_{LC}(\omega)&{\mathbf{G}}_{LR}(\omega)\\ {\mathbf{G}}_{CL}(\omega)&{\mathbf{G}}_{CC}(\omega)&{\mathbf{G}}_{CR}(\omega)\\ {\mathbf{G}}_{RL}(\omega)&{\mathbf{G}}_{RC}(\omega)&{\mathbf{G}}_{RR}(\omega)\end{array}\right]\;.$$

$${\mathbf{G}}_{CC}(\omega )=\frac{1}{(\omega +i\eta )\text{mathbbm}{1}_C-{\mathbf{H}}_{CC}-{\mathbf{\Sigma }}_L^{\mathrm{eq}}(\omega )-{\mathbf{\Sigma }}_R^{\mathrm{eq}}(\omega )},$$

$${\mathbf{\Sigma }}_{\alpha }^{\mathrm{eq}}(\omega )={\mathbf{H}}_{C\alpha }\frac{1}{(\omega +i\eta )\text{mathbbm}{1}_{\alpha }-{\mathbf{H}}_{\alpha \alpha }}{\mathbf{H}}_{\alpha C}.$$

$$\mbox{\boldmath$\rho$}^{\rm eq}=2\sum_{\alpha=L,R}\int\!\frac{{\rm d}\omega}{2\pi}f(\omega){\mathbf{G}}(\omega){\mathbf{\Gamma}}^{\rm eq}_{\alpha}(\omega){\mathbf{G}}^{\dagger}(\omega),$$

$${\mathbf{\Gamma }}_{\alpha }^{\mathrm{eq}}(\omega )=i\left({\mathbf{\Sigma }}_{\alpha }^{\mathrm{eq}}(\omega )-{\mathbf{\Sigma }}_{\alpha }^{\mathrm{eq},†}(\omega )\right).$$

$$\mbox{\boldmath$\rho$}=2\sum_{\alpha=L,R}\int\!\frac{{\rm d}\omega}{2\pi}f(\omega-V_{\alpha}){\mathbf{G}}(\omega){\mathbf{\Gamma}}_{\alpha}(\omega){\mathbf{G}}^{\dagger}(\omega)\;$$

$$n(\mathbf{r})=2\sum _{\alpha =L,R}\int \frac{\mathrm{d}\omega }{2\pi }f(\omega -V_{\alpha })A_{\alpha ,s}(\mathbf{r},\omega )$$

$$A_{\alpha ,s}(\mathbf{r},\omega )=\sum _{jq,j^{\prime }{q}^{\prime }}[{\mathbf{A}}_{\alpha ,s}(\omega ){]}_{jq,j^{\prime }{q}^{\prime }}⟨jq|\mathbf{r}⟩⟨\mathbf{r}|j^{\prime }{q}^{\prime }⟩,$$

$${\mathbf{A}}_{\alpha ,s}(\omega )\equiv \mathbf{G}(\omega ){\mathbf{\Gamma }}_{\alpha }(\omega ){\mathbf{G}}^{†}(\omega ).$$

$$N$$

$$I_s$$

$$G_s$$

$$\int {\mathrm{d}}^{3}rn(\mathbf{r})=N+\eta$$

$$E_{N+\eta }=(1-\eta )E_N+\eta E_{N+1}$$

$$I(N)=E_{N-1}-E_N,$$

$$A(N)=E_N-E_{N+1}.$$

$$\Delta (N)$$

$${\Delta }_s(N)$$

$$\Delta_{\rm xc}(N)=\lim_{\eta\to 0^{+}}\left(\frac{\delta E_{\rm xc}[n]}{\delta n({\bf r})}\bigg{|}_{N+\eta}-\frac{\delta E_{\rm xc}[n]}{\delta n({\bf r})}\bigg{|}_{N-\eta}\right)\;.$$

$$\Delta (N)={\Delta }_s(N)+{\Delta }_{\mathrm{xc}}(N).$$

$${\hat{H}}^{\mathrm{dot}}=v\hat{n}+U{\hat{n}}_{\uparrow }{\hat{n}}_{\downarrow }$$

Full text

I Introduction

Density Functional Theory (DFT) is probably the most popular method for an ab-initio description of atoms, molecules, and solids both in DreizlerGross:90 and out of (thermal) equilibrium Ullrich:12 ; Maitra:16 . The main reason for this popularity is its relative numerical simplicity which arises due to the mapping of the interacting many-electron problem onto an effective non-interacting one. This simplicity has given rise to a plethora of applications of DFT to the description of more and more complex systems in atomistic detail. However, there are also physical situations where DFT has been less successful, one of those being the description of strongly correlated systems. Since the fundamental theorems of DFT, both in thermal equilibrium and in its time-dependent out-of-equilibrium incarnation, give the framework a sound theoretical foundation, the failures of DFT have to be attributed to the inadequacy of known approximations to the famous exchange-correlation (xc) functional.

In this Topical Review we are concerned with a DFT treatment of the particular physical situation of quantum transport, i.e., the description of electronic transport through a nanoscale region contacted by metallic leads and driven out of equilibrium by application of an external bias. We will discuss different DFT frameworks for the description of quantum transport with a particular focus on recent advances in DFT approximations to deal with strongly correlated nanojunctions. It turns out that in all these frameworks, if one wants to capture strong correlations within DFT, the corresponding xc functionals need to have step features at integer occupation. These steps are intimately related to the famous derivative discontinuity of the exact xc potential.

The first application of DFT to electronic transport goes back to a seminal paper of Lang Lang:95 where, following ideas of Landauer Landauer:57 and Büttiker Buettiker:86 , transport in the steady state is treated as a scattering problem of effectively non-interacting electrons. The resulting formalism, known as Landauer-Büttiker plus DFT (LB+DFT), or its equivalent formulation in terms of non-equilibrium Green functions (NEGF) TGW.2001 ; TGW2.2001 ; BMOTS.2002 are by far the most widely used DFT schemes for transport. In Sec. II.1 we will present this formalism and critically discuss its merits and shortcomings. We will show that the LB+DFT framework, at zero temperature and in the limit of small bias, is capable of correctly capturing features related to the Kondo effect, a hallmark of strong electronic correlation. On the other hand, at finite temperatures LB+DFT fails to correctly describe Coulomb blockade, another ubiquitous correlation effect.

The Landauer-Büttiker formalism is only concerned with steady-state electronic transport under application of a DC bias voltage. Alternatively one may view transport as an explictly time-dependent problem where a system contacted to metallic leads is initially in thermal equilibrium. At a given time $t_{0}$, the system is driven out of equilibrium by application of a bias and one follows its time evolution C.1980 ; sa-1.2004 . In a DFT framework, this situations can naturally be described with Time-Dependent DFT (TDDFT) rg.1984 ; Ullrich:12 ; sa-2.2004 ; ksarg.2005 . Of course, application of TDDFT to transport is ideally suited to treat systems moving under the influence of explicitly time-dependent fields (e.g., AC bias). It turns out, however, that even for the steady state which develops in the long-time limit after application of a DC bias, TDDFT in principle leads to corrections to the LB+DFT formalism sa-2.2004 ; ewk.2004 ; kbe.2006 ; skgr.2007 . The TDDFT approach to transport will be discussed in Sec. III. The TDDFT corrections to LB+DFT are shown to be crucial in a correct description of Coulomb blockade in the zero-bias conductance. However, the price to be paid now is a lack of Kondo features in this framework.

Finally, in Sec. IV we review yet another, very recent DFT approach to steady-state transport called i-DFT StefanucciKurth:15 . Compared to the LB+DFT approach, the novelty is that the corresponding Kohn-Sham (KS) system is characterized not only in terms of the exchange-correlation (xc) potential in the nanojunction but also in terms of an xc contribution to the bias. For this framework we first construct approximations which are capable of correctly describing Coulomb blockade not only in the zero-bias limit but also at finite bias. Again, these approximations miss Kondo physics which, however, can be included using rather simple arguments. The i-DFT framework thus allows for a unified description of both Kondo effect and Coulomb blockade in finite-bias electronic transport as well as a transition between these regimes as temperature increases.

II Quantum Transport with Density Functional Theory

II.1 The standard approach: Landauer-Büttiker Formalism

In the typical setup of transport, a mesoscopic or nanoscopic central region $C$ is connected to two (or more) metallic leads and one is interested in the current flowing through $C$ upon application of a bias between any two leads.

For simplicity, here we only consider the case of the central region $C$ connected to left ($L$) and right ($R$) leads, see Fig. 1 for an illustration of the setup. A DFT description of such a system in thermal equilibrium at temperature $T=1/\beta$ and chemical potential $\mu$ requires the self-consistent solution of the Kohn-Sham (KS) equations (we use atomic units unless otherwise stated)

[TABLE]

where $v_{0}$ is the external potential generated by the positively charged nuclei and

[TABLE]

is the electronic density, $f(\epsilon)=1/(e^{\beta(\epsilon-\mu)}+1)$ being the Fermi function. In Eq. (1) there appear two functionals of the density. These are the Hartree potential

[TABLE]

with $w({\bf r}-{\bf r}^{\prime})$ the electron-electron interaction, and the exchange-correlation (xc) potential

[TABLE]

The xc potential is the functional derivative of the xc (free) energy functional $E_{\rm xc}[n]$.

A first technical difficulty when trying to solve the KS equations in the geometry of (semi-infinite) left and right leads connected to a central region arises from the fact that such a system is neither finite nor periodic. Thus standard techniques for electronic structure calculations of either finite or periodic systems need to be adapted to the problem.

One convenient way of dealing with the transport geometry is to use the language of Green’s functions. To begin with, we introduce a localized single-particle basis $|jq\rangle$ where $j$ denotes an atomic site and $q$ labels the different basis functions localized at this site. For simplicity, we assume that the single-particle basis is orthonormal and complete, i.e., $\langle jq|j^{\prime}q^{\prime}\rangle=\delta_{j,j^{\prime}}\delta_{q,q^{\prime}}$ and $\sum_{j,q}|jq\rangle\langle jq|=\mathbbm{1}$. In this basis, the KS Hamiltonian \mbox{\boldmath\mathcal{H}}_{s} can be written in $3\times 3$ block form as

[TABLE]

where ${\mathbf{H}}_{\alpha\alpha^{\prime}}$ collects all matrix elements connecting regions $\alpha$ and $\alpha^{\prime}$ ($\alpha,\alpha^{\prime}\in\{L,C,R\})$. Hereafter we use boldface letters to denote matrices in the one-electron basis. Note that in Eq. (5) we have assumed that all matrix elements connecting the left and right leads vanish. The KS single-particle orbitals $\psi_{k}$ can be expanded in the localized basis as

[TABLE]

which allows to write the density according to

[TABLE]

In Eq. (7) \mbox{\boldmath\rho}^{\rm eq} is the equilibrium one-particle density matrix whose matrix elements are given by

[TABLE]

where the factor of $2$ comes from spin.

The retarded Green’s function ${\mathbf{G}}(\omega)$ at energy $\omega$ is defined through

[TABLE]

with the infinitesimal $\eta\to 0^{+}$. We write the Green’s function in the same block structure as the Hamiltonian

[TABLE]

Using Eq. (9) we can easily solve for ${\mathbf{G}}_{CC}$ and find

[TABLE]

where the (retarded) embedding self energy for lead $\alpha$ is defined as

[TABLE]

Hereafter we will omit the subscript $CC$ from all matrices with both indices in the central region. Knowledge of ${\mathbf{G}}$ allows us to obtain the central block \mbox{\boldmath\rho}^{\rm eq} of the density matrix according to

[TABLE]

where the broadening matrix

[TABLE]

It is worth noting that ${\mathbf{G}}$ depends, through the Hartree-xc potential $v_{\rm Hxc}\equiv v_{\rm H}+v_{\rm xc}$, on the density and thus, via Eq. (7), on the density matrix. In principle, $v_{\rm Hxc}$ depends on the density both in the leads and the central region and therefore Eq. (13) is not a closed equation for \mbox{\boldmath\rho}^{\rm eq}. However, if any local or semilocal approximation such as LDA or GGA is employed and, at the same time, the embedding self energies ${\mathbf{\Sigma}}_{\alpha}$ for the leads are known, then Eq. (13) indeed becomes a self-consistent equation for \mbox{\boldmath\rho}^{\rm eq}.

So far, we have been discussing the situation in thermal equilibrium. In transport, however, one is interested in the scenario where the system is driven out of equilibrium by application of a bias. Most commonly, one is concerned with the steady state current of the system under application of a DC bias. In the picture of Landauer Landauer:57 and Büttiker Buettiker:86 , this steady current may be viewed as the result of lead electrons scattering off the potential of the central region $C$. In a seminal paper Lang:95 , Lang proposed to use the KS potential of DFT in the central region as the scattering potential. The idea is to calculate the scattering states deep in the left and right leads and filling them up to the chemical potentials $\mu_{\alpha}=\mu-V_{\alpha}$ ($\alpha=L,R$) shifted by the bias $V_{\alpha}$ in lead $\alpha$. One can now proceed by directly calculating the scattering states via, e.g., the Lippmann-Schwinger equation Lang:95 ; BMOTS.2002 . Equivalently, one may use the non-equilibrium Green’s function (NEGF) formalism svl-book to calculate the steady-state density matrix or Green’s function of the central region. The central block of the steady-state density matrix $\rho$ can then be obtained by the following equation

[TABLE]

which is structurally very similar to Eq. (13), the only difference being that ${\mathbf{G}}$ is calculated with \mbox{\boldmath\Sigma}_{\alpha}(\omega)=\mbox{\boldmath\Sigma}_{\alpha}^{\rm eq}(\omega-V_{\alpha}), the broadening matrix \mbox{\boldmath\Gamma}_{\alpha}(\omega)=\mbox{\boldmath\Gamma}_{\alpha}^{\rm eq}(\omega-V_{\alpha}) and the Fermi function for lead $\alpha$ is shifted by $V_{\alpha}$. With the density matrix one can then calculate the electronic density as in Eq. (7) using $\rho$ in place of \mbox{\boldmath\rho}^{\rm eq}. After some elementary algebra one finds

[TABLE]

where

[TABLE]

and

[TABLE]

is the partial KS spectral function. Again, through the dependence of the Green’s function on the Hxc potential and thus on the density, this defines a self-consistency problem. From the local density we can also calculate the total number of electrons in region $C$ as

[TABLE]

where we have taken into account that the states $jq$ form a complete set in region $C$ and the trace is over the single-particle basis in region $C$ only.

The steady current can be calculated using the self-consistent Green’s function via the celebrated Landauer-Büttiker (LB) formula

[TABLE]

The subscript $s$ in $I_{s}$ highlights the fact that Eq. (20) gives the steady current of the KS system. There exist no formal proof that $I_{s}$ is the same as the steady current $I$ of the interacting system. In fact, in this review we will present relevant cases for which $I_{s}\neq I$. From Eq. (20) we can also calculate the KS zero-bias conductance

[TABLE]

where all quantities in the trace are evaluated at $V_{\alpha}=0$, i.e., at equilibrium.

The formalism described above combines the LB formula with DFT and it is widely used to describe steady-state transport through nanoscale systems with atomistic detail. This level of description is essential in the field of molecular electronics me-book1 ; me-book2 whose central tenet is to use single molecules as active electronic devices.

II.2 The derivative discontinuity of the DFT exchange-correlation potential

In typical applications of DFT one aims to calculate the ground state energy and/or the electronic structure of a molecule or solid for a given, fixed number of electrons. However, in the transport problem a nanostructure such as a quantum dot or a molecule is connected to leads and the number of electrons in the nanostructure fluctuates. In other words, instead of dealing with a closed system at fixed particle number, in transport we are studying open systems connected to particle reservoirs. Therefore one may expect that the dependence of density functional approximations on the particle number becomes important. This will be the concern of the present section.

In a seminal paper PerdewParrLevyBalduz:82 , Perdew and coworkers pointed out that the exact xc potential of DFT at zero temperature exhibits discontinuous steps as the particle number crosses an integer. To show this, they constructed an ensemble DFT based on ensembles of states with different (integer) electron numbers. The ensemble expectation value of the electron number operator then can yield non-integer values, or, in other words, the (ensemble) density integrates to a non-integer

[TABLE]

with $N\in\mathbbm{N}$ and $0\leq\eta<1$. Perdew and coworkers proved that the ground ensemble energy is a piecewise linear function of the (fractional) electron number, i.e.,

[TABLE]

where $E_{N}$ is the ground state energy at integer particle number $N$. Therefore, the ground ensemble energy as function of electron number is given by a series of straight lines connecting the ground state energies with integer numbers of particles. The slope of these straight line segments may be expressed in terms of physical quantities: the ionization potential $I(N)$ of the $N$-electron system is

[TABLE]

while the electron affinity $A(N)$ is given by

[TABLE]

As a consequence, the discontinuous jump of the derivative of the ground ensemble energy at integer electron number can be expressed as

[TABLE]

Decomposing the total energy functional into its components, there are only two terms which are discontinuous as the particle number crosses an integer. The first term comes from the non-interacting kinetic energy $T_{s}[n]$, i.e.,

[TABLE]

where $\varepsilon_{k}(N)$ is the $k$-th lowest KS energy eigenvalue of an $N$ electron calculation. The second terms is the famous derivative discontinuity of the xc potential

[TABLE]

The total discontinuity is the sum of the two, i.e.,

[TABLE]

The xc discontinuity not only gives an important contribution to the fundamental gap of semiconductors and insulators, but it is also a highly relevant property of the exact xc energy functional in other situations. For instance, it is exactly the xc discontinuity which ensures that heteronuclear molecules dissociate into fragments with integer electron numbers Perdew:85-2 ; RuzsinszkyPerdewCsonkaVydrovScuseria:06 ; FuksMaitra:14 . The development of the xc discontinuity from solvable systems with fractional electron number as the fraction $\eta$ approaches zero has been studied in Refs. SagvoldenPerdew:08 ; GoriGiorgiSavin:09 . The importance of the derivative discontinuity has emerged in other contexts too. For instance, the exact solution of the one-dimensional Hubbard model shows that at zero temperature the (uniform) xc potential is discontinuous as the number of particles per site crosses unity (half-filling) LimaSilvaOliveiraCapelle:03 ; KVOC.2011 ; CapelleCampo:13 . In more than one-dimension it has been shown that the discontinuity occurs only for sufficiently strong interactions and that the critical value of the interaction strength is the same as that of the Mott-Hubbard transition KarlssonPriviteiraVerdozzi:11 ; KKPV.2013 .

Popular density functional approximations such as LDA or GGA’s exhibit a vanishing derivative discontinuity, at least in the way they are commonly used. It has been pointed out only recently, however, that even for these functionals one can construct derivative discontinuities KraislerKronik:13 : to do so one interprets the corresponding energy functionals as orbital functionals in an ensemble DFT framework and consequently calculates the xc potential within the optimized effective potential framework. There are also density functional approximations which do exhibit a derivative discontinuity, most prominently perhaps explicitly orbital dependent functionals such as the exact-exchange functional. Other density functional approximations based on best-fitting are available for the xc potential of the one-dimensional Hubbard model LimaSilvaOliveiraCapelle:03 ; CapelleCampo:13 . In this case one can also show how the discontinuous steps emerge in the low temperature limit xctk.2012 .

In the next Section we consider a different class of interacting models showing a discontinuity in the xc potential. Then, in Sec. II.4, we will use the LB+DFT approach to study this class of models in a quantum transport setup.

II.3 The derivative discontinuity in a few illustrative examples

The importance of the derivative discontinuity appears in the description of systems whose electron number can fluctuate. This is precisely the situation of junctions (the central region $C$ of Sec. II.1) connected to leads. In the limit of vanishing contacts the equilibrium properties of the junction are given by the grand canonical partition function. Below we will use the grand canonical generalization of DFT by Mermin m.1965 in two paradigmatic model systems and provide a somewhat different perspective on the derivative discontinuity.

II.3.1 Single Site Hubbard Model

Our first model for a quantum dot consists of a single level with on-site energy $v$ which can hold up to two electrons sk.2011 . It is described by the Hamiltonian

[TABLE]

where $U$ is the charging energy, $\hat{n}_{\sigma}=\hat{d}^{\dagger}_{\sigma}\hat{d}_{\sigma}$ is the number operator for electrons of spin $\sigma$ on the dot, $\hat{n}=\hat{n}_{\uparrow}+\hat{n}_{\downarrow}$, and $\hat{d}^{\dagger}_{\sigma}$ and $\hat{d}_{\sigma}$ are the corresponding electron creation and annihilation operators, respectively.

The KS Hamiltonian of such a system is

[TABLE]

According to Mermin’s finite-temperature version of DFT m.1965 there exists a unique potential $v_{s}$ for which this KS Hamiltonian yields, for a given temperature and chemical potential, the same density as the interacting Hamiltonian. We write $v_{s}=v+v_{\rm Hxc}[n]$ where the Hxc potential $v_{\rm Hxc}[n]$ of the single site model (SSM) is a function of the dot density $n$. Let us derive the exact functional form of the SSM Hxc potential.

We start by observing that for both the interacting and KS Hamiltonian, see Eqs. (30) and (31), the eigenstates for electron occupation zero, one, or two are: $|0\rangle$, $|\uparrow\rangle$, $|\downarrow\rangle$, and $|\uparrow\downarrow\rangle$. The eigenvalues of $\hat{H}^{\rm dot}$ corresponding to these states are [math], $v$, $v$, and $2v+U$, respectively, while for $\hat{H}^{\rm dot}_{s}$ the eigenvalues are [math], $v_{s}$, $v_{s}$, and $2v_{s}$. We now consider the single site at thermal equilibrium in contact with a bath at inverse temperature $\beta=1/T$ and chemical potential $\mu$. The total density on the interacting dot reads

[TABLE]

where we have defined $\tilde{v}=v-\mu$ and the partition function

[TABLE]

For the KS dot the result simply is

[TABLE]

with $\tilde{v}_{s}=v_{s}-\mu$. Both density-potential relations (32) and (34) can be inverted analytically which allows us to write the exact Hxc potential of the single-site model as

[TABLE]

where

[TABLE]

This is an odd function of its argument, $g_{U}(-x)=-g_{U}(x)$, and therefore $v_{\rm Hxc}[n=1]=\frac{U}{2}$ for all temperatures. In the left panel of Fig. 2 we show the Hxc potential (35) for various temperature. The most prominent feature is the (smoothened) step at half-filling ($n=1$) and low temperatures. In the zero temperature limit, the Hxc potential approaches a discontinuous step function $v_{\rm Hxc}(n)\stackrel{{\scriptstyle T\rightarrow 0}}{{\longrightarrow}}U\theta(n-1)$ where $\theta(x)$ is the Heaviside step function. This step is nothing but the derivative discontinuity at zero temperature discussed above. Hence, the discontinuity naturally emerges in the zero-temperature limit of Mermin’s formulation of DFT. At finite temperature, on the other hand, the discontinuity is thermally broadened. Thermal broadening is not the only mechanism to smoothen the step in the xc potential. As we will discuss below, particularly relevant for transport is the broadening due to the presence of contacts.

As a pedagogical illustration of the consequences of the step feature in the xc potential, we calculate the density of the single site as function of the on-site potential $v$ from solving the exact KS self-consistency condition

[TABLE]

as well as the KS self-consistency condition at the Hartree level

[TABLE]

where the SSM Hartree potential

[TABLE]

is a linear (and hence continuous) function of $n$. The right panel of Fig. 2 shows the results (for convenience we choose $\mu=0$). The most striking difference is the plateau feature of the exact density in the low-temperature regime for gate potentials in the range $-U\lesssim v\lesssim 0$. This is completely missing in the Hartree approximation (or any LDA or GGA for that matter) since it is a direct consequence of the step in $v_{\rm Hxc}$. For increasing temperatures, on the other hand, the Hartree density approaches the exact one as the step feature is washed out.

II.3.2 Constant Interaction Model

Our second model is slightly more complicated and is known as the Constant Interaction Model (CIM). The CIM Hamiltonian reads

[TABLE]

where $\epsilon_{i}$ are the energies of the single-particle levels and $U$ is the repulsive energy between two electrons. For $M=1$ and $\varepsilon_{1}=v$, the CIM Hamiltonian reduces to the SSM Hamiltonian of Eq. (30). Notice that the interacting part of the CIM Hamiltonian can be equivalently written as

[TABLE]

where $\hat{N}=\sum_{i\sigma}\hat{n}_{i\sigma}$ is the operator of the total number of particles.

The KS Hamiltonian of the CIM is

[TABLE]

with the Hxc potential $v_{\rm Hxc}$ and the density operator $\hat{n}({\bf r})$ which is related to $\hat{N}$ by $\hat{N}=\int\!\!{\rm d}^{3}r\;\hat{n}({\bf r})$. It can be shown StefanucciKurth:13 that in the limit of zero temperature, the exact Hxc potential of the CIM is independent of position and depends only on the total number $N$ of electrons in the system. In other words

[TABLE]

Let us prove this statement. We denote by $|\Psi^{N}_{l}\rangle$ the linearly independent ground states with $N$ electrons, where the index $l=1,\ldots,d_{N}$ and $d_{N}$ is the degeneracy of the ground-state multiplet of energy $E(N)$. In the zero-temperature limit of the grand canonical ensemble at chemical potential $\mu$, the number of electrons $N$ is the largest integer for which the addition energy fulfills

[TABLE]

The corresponding ground state density reads

[TABLE]

Due to the particular form of the interaction $\hat{H}_{\rm int}$, the eigenstates of the CIM Hamiltonian are the same as the eigenstates of $\hat{H}_{0}$; hence they are many-body Slater determinants with every level occupied by either zero or one electron of spin $\sigma$. This implies that

[TABLE]

where $n_{i\sigma}(N,l)=1$ if $|\Psi^{N}_{l}\rangle$ contains an electron of spin $\sigma$ on level $i$ and zero otherwise. We turn now to the KS Hamiltonian in Eq. (42). If the Hxc potential is independent of position then the CIM Hamiltonian and the KS Hamiltonian have the same eigenstates with the same degeneracies. Therefore the ground state densities in the $N$-electron sector are the same both in the interacting and in the KS system. Consequently, for given chemical potential $\mu$, the only role of the Hxc potential $v_{\rm Hxc}$ is to ensure that at zero temperature both the interacting and the KS system have the same number of electrons. This is achieved if $N$ is the largest integer such that

[TABLE]

Thus, for any real $N$, the explicit form of the CIM Hxc potential at zero temperature is

[TABLE]

where $\bar{N}={\rm Int}[N]$ is the integer part of $N$ and we have used Eq. (44) in the last step. We conclude that the CIM Hxc potential is piecewise constant with discontinuities of height $U$ whenever $N$ crosses an integer.

At finite temperature, the CIM Hxc potential in general is more complicated. However, if all the single-particle levels $\varepsilon_{i}$ are degenerate, it is again possible to write it as a position-independent constant depending only on the total number $N$ in the system. Note that now $N={\rm Tr}\left\{\exp(-\beta(\hat{H}^{\rm CIM}-\mu))\hat{N}\right\}/{\rm Tr}\left\{\exp(-\beta(\hat{H}^{\rm CIM}-\mu))\right\}$ is to be understood as the thermal average of the particle number operator $\hat{N}$.

We have constructed the finite temperature Hxc potential of a 6-level CIM with degenerate single-particle levels by numerical reverse engineering where we have used the scheme described in A to calculate the equilibrium occupation of the CIM. The resulting Hxc potential of this construction is shown in Fig. 3. The step structure at low temperature is apparent and in the zero-temperature limit approaches our result of Eq. (48), as it should be. As in the case of the isolated site, at high temperatures the steps are washed out and the Hxc potential approaches a linear function of $N$.

II.4 Successes and limitations of the LB+DFT formalism for strongly correlated transport

II.4.1 The single impurity Anderson model

In Sec. II.3.1 we have considered an uncontacted single-level quantum dot. In order to study electronic transport through this dot we connect it to a left ($L$) and right ($R$) lead described as one-dimensional, semi-infinite tight binding chains. The resulting model is the celebrated single-impurity Anderson model (SIAM) Anderson:61 for which the Hamiltonian reads

[TABLE]

Here, the tight-binding lead $\alpha=L,R$ is described by

[TABLE]

and the tunneling Hamiltonian connecting the dot to the leads is

[TABLE]

where we have assumed symmetric coupling to left and right leads. The embedding self energy for one-dimensional tight-binding leads is known analytically. However, here we concentrate on half-filled leads in the parameter regime $V_{\rm link}\ll V$, the so-called wide-band limit (WBL). In this limit the embedding self energy for lead $\alpha$ becomes a purely imaginary, energy independent constant $\gamma_{\alpha}=2V_{\rm link}^{2}/V$ and the only relevant energy scale for electron tunneling is $\gamma=\gamma_{L}+\gamma_{R}$.

For the SIAM the single-particle density matrix $\rho$ is a $1\times 1$ matrix whose only entry is the value of the density on the dot. Therefore, the self-consistency condition of Eq. (15) becomes a nonlinear equation for the only unknown \mbox{\boldmath\rho}=n. Taking advantage of the WBL nature of the leads it is easy to show that in thermal equilibrium (no bias) Eq. (16) becomes

[TABLE]

This equation correctly reduces to Eq. (37) in the limit of vanishing contacts, i.e., $\gamma\to 0$. While for the isolated dot the analytic form of the xc potential is known exactly, see Eq. (35), this is not the case for the contacted dot. However, one may still use $v_{\rm Hxc}[n]$ of Eq. (35) as an approximation to the exact xc potential if $\gamma\ll T$. With the self-consistent density $n$ solving Eq. (52) we can easily compute the KS zero-bias conductance via Eq. (21)

[TABLE]

where $G_{0}=1/\pi$ is the quantum of conductance.

In Fig. 4 we show $G_{s}$ for two different temperatures and compare it to the SIAM conductance $G$ as obtained from numerical renormalization group (NRG) calculations IzumidaSakaiSuzuki:01 ; IzumidaSakai:05 . In the left panel, the temperature is much smaller than the Kondo temperature JakobsPletyukhovSchoeller:10

[TABLE]

whereas in the right panel $T\gg T_{\rm K}$. Below the Kondo temperature the most prominent feature of $G_{s}$ and $G$ is the plateau at one quantum of conductance $G_{0}$ in the region of gate potentials $-U\lesssim v\lesssim 0$. Physically this plateau is due to the screening of the single electron spin at the impurity by a spin cloud of opposite-spin electrons at the interface of both leads, a phenomenon known as Kondo effect. In the context of DFT, the plateau originates from the pinning of the the KS level to the Fermi energy which is a direct consequence of the use of an Hxc potential with a step at half filling. This is confirmed by the absence of the conductance plateau in DFT calculations using the Hartree potential of Eq. (39). Furthermore, we observe that in the DFT calculation using $v_{\rm Hxc}$, the conductance plateau extends further and terminates more abruptly than in the reference NRG calculation.

For $T\gg T_{\rm K}$, see right panel of Fig. 4, DFT even qualitatively disagrees with the NRG result. While the NRG conductance now clearly shows two peaks at $v\approx-U$ and $v\approx 0$ due to Coulomb blockade, the DFT calculation still exhibits a conductance plateau, although at a value smaller than $G_{0}$. Below we discuss in some detail the reasons for the success of the LB+DFT approach at low temperatures and its failure at high ones.

At first sight it may come as a surprise that features of complicated many-body physics such as the Kondo effect can qualitatively be captured with a simple DFT model and an explanation is called for. A first explanation can be gleaned from the Meir-Wingreen formula of the zero bias conductance MeirWingreen:92

[TABLE]

with $A(\omega)$ the interacting spectral function. At $T=0$ the Meir-Wingreen formula gives

[TABLE]

where ${\mathbf{G}}^{-1}(\omega)=[\omega-v-\Sigma(\omega)+i\frac{\gamma}{2}]$ is the $1\times 1$ interacting Green’s function and $\Sigma(\omega)$ is the many-body self-energy. At the Fermi energy, quasiparticles have an infinite lifetime, i.e., we have ${\rm Im}\,\Sigma(\mu)=0$. Therefore one can see from Eq. (56) that it is indeed possible to reproduce the exact conductance from a KS calculation if the KS potential at the impurity is $v_{s}=v+{\rm Re}\,\Sigma(\mu)$. A second explanation that LB+DFT can give the exact zero-bias conductance can be found in the Friedel sum rule Langreth:66 ; mera-1 ; mera-2 which states that at zero temperature the zero-bias conductance of the SIAM is fully determined by the ground state density at the impurity. Since exact DFT by construction gives the exact ground state density, it therefore also must yield the exact zero-bias conductance, including the conductance plateau due to the Kondo effect. The argument of the Friedel sum rule in the context of the LB+DFT conductance for the SIAM has independently been discussed in Refs. sk.2011 ; blbs.2012 ; tse.2012 .

The above arguments also make clear where the small difference between the NRG and KS conductances at $T\ll T_{\rm K}$ come from. The Hxc potential used in our calculation does not contain any information related to the contact to the leads. In other words the $v_{\rm Hxc}$ of Eq. (35) is the exact SIAM Hxc potential only for $\gamma\to 0$. However, the contacts are responsible for a broadening of the step structure and hence $v_{\rm Hxc}$ should exhibit a smeared step for finite $\gamma$ even at $T=0$. To include the broadening due to the contacts in the Hxc potential we proceed as follows. The exact spectral function of the isolated dot reads

[TABLE]

Broadening the delta-function $\delta(\omega)$ into a Lorentzian $\ell_{\gamma}(\omega)=\gamma/(\omega^{2}+\gamma^{2}/4)$ of width $\gamma$ we obtain the model many-body (MB) spectral function

[TABLE]

With this model spectral function the interacting density as function of $\tilde{v}=v-\mu$ is

[TABLE]

which can be inverted numerically to give $\tilde{v}[n]$. Similarly, for a non-interacting impurity with on-site potential $v_{s}$, through Eq. (59) for $U=0$ the density becomes a function of $\tilde{v}_{s}=v_{s}-\mu$ which again can be inverted numerically to give $\tilde{v}_{s}[n]$. From these results we then obtain a numerical model xc potential for the coupled impurity through

[TABLE]

This Hxc potential reduces to the $v_{\rm Hxc}$ of Eq. (35) for $\gamma\to 0$ and it has the nice feature of smearing the step at half filling with a width $\approx\gamma$ for $T\ll\gamma$ and a width $\approx T$ for $T\gg\gamma$. The quantitative effects of thermal broadening versus contact broadening are illustrated in Fig. 5. We wish to emphasize that independently of the nature of the broadening the qualitative behavior of the density does not change: in both cases $n$ exhibits a plateau of height $1$ (half-filling) as function of $v$ since $v_{\rm Hxc}$ pins the KS level to the chemical potential.

At zero temperature a reasonable fit to the model Hxc potential of Eq. (60) is given by the expression

[TABLE]

where the fitting parameter $W$ is defined as

[TABLE]

We will use this form at various places throughout this review. The (numerically) exact Hxc potential for the Anderson model at zero temperature has been constructed by a reverse engineering procedure using DMRG methods in Refs. blbs.2012, ; LiuBergfieldBurkeStafford:12, where the authors also give accurate parametrizations of their results.

We now turn to the discussion of our results for temperatures $T\gg T_{\rm K}$. From the right panel of Fig. 4 we see that the KS conductance even qualitativeley disagrees with the correct NRG result: instead of the Coulomb blockade peaks the DFT results still show a plateau, although at values smaller than $G_{0}$. In fact, both arguments given to explain the correct DFT description at zero temperature fail at finite temperature. When calculating the conductance from the Meir-Wingreen formula at finite temperature, we cannot restrict the discussion of the many-body Green’s function to the Fermi energy alone and thus the argument given above does not apply. Similarly, the validity of the Friedel sum rule is restricted to zero temperature. Therefore it is not surprising that the DFT conductance does not give the correct physics. This can explicitly be confirmed by considering a special value for the gate potential, the so-called particle-hole (ph) symmetric point at $v=-U/2$. At this value of the gate, the impurity is at half-filling ($n=1$) for all temperatures. In the DFT framework this means that the exact Hxc potential must have the value $v_{\rm Hxc}[n=1]=\frac{U}{2}$ such that the total KS potential vanishes. Our model Hxc potential in Eq. (60) correctly satisfies this condition, i.e., the value of $v_{\rm Hxc}$ for $n=1$ is exact for all temperatures and interaction strengths $U$. At the ph symmetric point, the conductance $G^{\rm ph}$ of the interacting Anderson impurity is a universal function of $T/T_{\rm K}$ which is known numerically Costi:00 ; AleinerBrouwerGlazman:02 . Therefore at the ph symmetric point we can compare the exact KS conductance $G^{\rm ph}_{s}$ with the exact $G^{\rm ph}$. The result is shown in Fig. 6. We see that while the exact KS conductance is correct at zero (and very low) temperature, it is widely off the mark over a wide temperature range.

II.4.2 Constant Interaction Model coupled to leads

In the previous Section we have seen what the LB+DFT formalism can and cannot describe for the Anderson model. In the present one we move to multi-level systems, in particular systems described by the CIM of Sec. II.3.2. There we have seen that the Hxc potential at zero temperature is independent of position and depends only on the total number $N$ of electrons on the (multi-level) quantum dot. Thus, in the basis of the single-particle eigenstates, the CIM KS Hamiltonian of Eq. (42) can be written as

[TABLE]

where $\hat{n}_{i}=\sum_{\sigma=\uparrow,\downarrow}\hat{n}_{i\sigma}$ and $M$ is the number of levels. In order to study the conductance properties of this system we have to connect it to leads. The non-interacting leads are described by the Hamiltonian

[TABLE]

where the $\hat{c}_{k\sigma\alpha}$ ($\hat{c}^{\dagger}_{k\sigma\alpha}$) create (annihilate) an electron with energy $\varepsilon_{k\alpha}$ and spin $\sigma$ in lead $\alpha$. The tunneling Hamiltonian takes the usual form, i.e.,

[TABLE]

which, from Eq. (14), leads to the broadening matrix $\Gamma_{ii^{\prime}}=2\pi\sum_{k\sigma\alpha}T_{k\alpha,i}^{*}T_{k\alpha,i^{\prime}}\delta(\omega-\varepsilon_{k\alpha})$. We will again work in the WBL where the $\Gamma_{ii^{\prime}}$ become independent of frequency. We further consider the $\Gamma_{ii^{\prime}}$ much smaller than the level spacings so that we can approximate

[TABLE]

in the dot. This assumption is certainly justified at low enough temperatures. At low temperatures we also expect that the essential qualitative feature of $v_{\rm Hxc}$ is a series of smeared steps at integer $N$. On this ground, the approximation we propose consists in summing the SIAM Hxc potential over all possible charged states of the CIM  PS.2012 . Taking into account Eq. (61) the explicit form of the approximation reads

[TABLE]

where, for later use, we already allowed for charging energies $U_{J}$ and widths $W_{J}$ which depend on the charging state $J$. In the present Section we consider $\gamma_{i}=\gamma$ independent of $i$ and use both charging energies $U$ and level broadenings $W=0.16\;\gamma/U$ (see Eq. (62)) independent of the charging state $J$.

We study the zero-bias conductance through the multi-level CIM. It is therefore sufficient to solve the KS problem in equilibrium. Taking into account that the Hxc potential depends only on the total $N$, there is only one self-consistency equation

[TABLE]

with \mbox{\boldmath\Gamma}=\mbox{\boldmath\Gamma}_{L}+\mbox{\boldmath\Gamma}_{R} and the trace is over all single-particle states of the quantum dot, see Eqs. (19). According to Eq. (11) the KS Green’s function is given by

[TABLE]

where the KS single-particle Hamiltonian has matrix elements $[{\mathbf{H}}]_{jj^{\prime}}=\delta_{jj^{\prime}}(\varepsilon_{j}+v_{\rm Hxc}[N])=:\delta_{jj^{\prime}}\varepsilon_{s,j}$ with the KS single-particle energies $\varepsilon_{s,j}$. At self-consistency we can calculate the KS zero-bias conductance using Eq. (21).

As a first example of transport through multiple correlated levels described by the CIM we consider two spin-degenerate single-particle levels, the HOMO and the LUMO, coupled to two wide-band leads. As functions of gate potential $v$, the KS HOMO and LUMO leves are given by

[TABLE]

and

[TABLE]

In Fig. 7 we show zero temperature results for density, KS conductance, and KS single-particle levels for the following parameters: HOMO-LUMO splitting $\Delta\varepsilon=0.5$ and $\gamma=0.05$ (all energies are given in units of $U$). At large positive gate potentials, the quantum dot is unoccupied (see upper panel). As the gate is lowered, the first electron enters the dot when the (KS) HOMO level becomes lower than the Fermi energy $\varepsilon_{\rm F}$. Due to the interaction, the second electron is blocked from entering the dot at the same gate and only can enter when the gate is lowered by $U$. In order for the third electron to enter the dot, we have to lower the gate by an additional energy $\Delta\varepsilon+U$, i.e., by the single-particle energy difference of the levels plus one additional charging energy. Finally, the fourth electron can only enter once the gate is lowered by another charging energy $U$. The KS conductance (middle panel) at zero temperature shows two regions of gate potentials with one quantum of conductance $G=G_{0}$. These regions correspond to those values of $v$ with an odd number of electrons occupying the quantum dot while for an even number of electrons the conductance is blocked, exactly the correct behavior for the appearance of the Kondo effect in multi-level quantum dots BruusFlensberg:04 . The way how this is achieved within the LB+DFT formalism can be deduced by looking at the KS energy levels (lower panel). We see that for those gate potentials for which the KS conductance equals $G_{0}$, one of the two KS levels is pinned to the Fermi energy $\varepsilon_{\rm F}$. Just like in the SIAM, this pinning leads to one of the conductance channels being open, i.e., $G=G_{0}$, and, of course, is a direct consequence of the step feature in the Hxc potential at integer occupation. At occupation with even number of electrons, none of the KS levels pins to $\varepsilon_{\rm F}$ and the conductance is blocked, despite the fact that $v_{\rm Hxc}$ also has steps at even integers of the occupation. As expected from our experience with the SIAM, at finite temperature the KS conductance qualitatively keeps its $T=0$ shape but with the plateau values now lowered. It completely fails to describe the transition from the Kondo to the Coulomb blockade regime.

As our second example we again consider the HOMO-LUMO model contacted to two leads but now for the case of degenerate HOMO and LUMO levels ($\Delta\varepsilon=0$) and all other parameters as before. The results are shown in Fig. 8. Again, the occupation exhibits plateaus (as function of gate) with integer occupation of the dot. The LB+DFT conductance (middle panel), however, now has a different structure. While for odd numbers of electrons on the dot, the conductance is still $G_{0}$, for occupation of two we have $G=2G_{0}$. This can be easily understood since in this case the whole dot is half filled, i.e., both (degenerate) single-particle levels are pinned to $\varepsilon_{\rm F}$ (see lower panel) and thus two conductance channels are open simultaneously. For those values of $v$ where $G=G_{0}$, the KS level doesn’t pin to $\varepsilon_{\rm F}$ but to $\pm\gamma/2$ (note the scale on the $y$-axis of the lower panel).

The examples discussed in this Sec. demonstrate that the LB+DFT approach correctly predicts the zero-temperature Kondo plateaus in $G$ as well as the height of the plateaus for degenerate levels. Of course, the success of LB+DFT strongly relies on the use of accurate Hxc potentials whose most important feature is the occurrence of steps as the number of particles crosses an integer. However, it is also clear that, just as in the SIAM, by increasing the temperature the transition from the Kondo to the Coulomb blockade regime is beyond the capabilities of the LB+DFT approach StefanucciKurth:13 . How to go beyond this approach will be the topic of the next two sections. Before closing this Section we wish to observe that in Figs. 7 and 8 we used the Hxc potential of Eq. (67) even for $T>\gamma$. As we discussed below Eq. (60), thermal broadening dominates for $T>\gamma$ and it would be more appropriate to use the SSM Hxc potential of Eq. (35). This is precisely what we did in the SIAM and, as we can see from Fig. 4, the only difference is that the width of the flanks of the Kondo plateau is proportional to $T$ instead of $\gamma$. Thus, the qualitative behavior of the total $N$ and of the KS conductance is independent of the nature (thermal or due to the contacts) of the broadening.

III Quantum Transport with Time-Dependent Density Functional Theory

The LB+DFT approach combines the noninteracting steady-state formulation of Landauer and Büttiker with the ground-state DFT formulation of interacting systems. This empirical combination suffers from a conceptual problem, i.e., the use of an equilibrium Hxc potential in a nonequilibrium steady-state situation. One way to solve this conceptual problem while still remaining in a one-particle framework is to reformulate the theory of quantum transport using the Time Dependent (TD) version of DFT. TDDFT allows for studying the TD current and density generated by an arbitrary TD bias and hence, as a special case, the real-time evolution of the current and density after the switch-on of an external DC bias. Steady-state quantities can simply be obtained as the long-time limit of the TDDFT results.

We consider again the geometry of the previous Section with a central region connected to left and right electrodes. The whole system is initially, say at time $t\leq 0$, in equilibrium at a certain temperature and chemical potential. The charge density of the electrodes is perfectly balanced and no current flows through the junction. The system is driven out of equilibrium by an external electric field, like the one generated by a battery. The formation of dipole layers screens the potential drop along the electrodes and the total potential turns out to be uniform in the left and right bulks. Accordingly, the potential drop, or bias, is confined in the neighborhood of the junction.

Let $\psi_{k}({\mathbf{r}})$ be a KS eigenstate of the equilibrium system with KS energy $\varepsilon_{k}$. The time evolution of this state is governed by the time-dependent KS equation

[TABLE]

where $v_{\rm ext}$ is the external potential, i.e., the sum of the nuclear potential and the time-dependent potential of the battery, and $v_{\rm Hxc}=v_{\rm H}+v_{\rm xc}$ is the Hartree-xc potential of TDDFT. The exact Hxc potential $v_{\rm Hxc}({\mathbf{r}},t)$ in point ${\mathbf{r}}$ at time $t$ depends on the density

[TABLE]

in all points ${\mathbf{r}}^{\prime}$ and at all times $t^{\prime}<t$. Hence Eq. (72) for all single-particle indices $k$ together with Eq. (73) form a nonlinear system of coupled differential equations.

The solution of the differential equations is, in general, a difficult task since (i) the system is infinitely large (due to the electrodes) and spatially non-periodic (due to the junction) and (ii) the KS wavefunctions at time $t=0$ are delocalized all over the system. Nevertheless, it is still possible to develop efficient numerical algorithms if we take advantage of the uniformity of the time-dependent part of the total potential in the left and right bulks. The property of uniformity means that

[TABLE]

where $V_{\alpha}(t)$ is the (experimentally measured) bias and

[TABLE]

where $V_{\alpha,\rm xc}(t)$ is the xc bias correction predicted by TDDFT. In Eqs. (74) and (75) the limit ${\mathbf{r}}\to\alpha$ signifies that ${\mathbf{r}}$ should be taken deep inside lead $\alpha$. Thus, the time-dependent part of the KS Hamiltonian in Eq. (72) is spatially constant in the leads.

In the last decade several algorithms have been proposed to obtain the one-particle density matrix $\rho({\mathbf{r}},{\mathbf{r}}^{\prime},t)=\sum_{k}f(\varepsilon_{k})\psi^{\ast}_{k}({\mathbf{r}},t)\psi_{k}({\mathbf{r}}^{\prime},t)$ from which to extract the density $n({\mathbf{r}},t)=\rho({\mathbf{r}},{\mathbf{r}},t)$, see Eq. (73), and the longitudinal current I(t)=\int d{\mathbf{r}}_{\perp}\,{\rm Im}\,[(\mbox{\boldmath\nabla}-\mbox{\boldmath\nabla}^{\prime})\rho({\mathbf{r}},{\mathbf{r}}^{\prime},t)]_{{\mathbf{r}}={\mathbf{r}}^{\prime}} (the integral is over a surface perpendicular to the current flow). Among the methods based on the solution of Eq. (72) we mention the wavefunction approach with transparent boundaries ksarg.2005 ; vsa.2006 ; spc.2008 ; spc.2010 ; kskvg.2010 and with absorbing boundaries bsin.2004 ; v.2011 ; gwshgw.2014 ; ww.2016 ; sk.2016 , the microcanonical approach bsdv.2005 ; cevv.2006 ; psc.2008 ; psc.2009 , the supercell approach qlly.2006 and the stroboscopic approach bcg.2008 ; kb.2014 . Other methods are based on the solution of the Dyson equation for the Green’s function mwg.2006 ; zwymc.2007 ; yzcwfn.2011 ; mgm.2007 ; pd.2008 ; cs.2009 ; PopescuCroy:16 ; zetal.2010 ; zcw.2013 ; tpsvl.2014 ; cuansing:17 , on the propagation of the Green’s function through the Kadanoff-Baym equations mssvl.2008 ; mssvl.2009 ; pva.2009 ; pva.2010 , on the calculation of Bohm trajectories o.2007 and on different types of master equations for the density matrix sssbht.2006 ; ly.2007 ; xetal.2012 ; zkh.2014 ; chf.2014 .

In Refs. sa-1.2004 ; sa-2.2004 it has been shown that if the bias $V_{\alpha}$ is constant for large times and if the KS system attains a steady state in the long time limit then the steady current flowing through the junction is given by

[TABLE]

The main difference between Eq. (76) and the LB+DFT formula in Eq. (20) is the appearence of the KS bias $V_{\alpha,s}$ in the Fermi functions. The KS bias is defined according to

[TABLE]

and differs from the physical bias by the addition of an xc correction. Another important difference between Eq. (76) and the LB+DFT formula in Eq. (20) lies in the calculation of the Green’s function ${\mathbf{G}}$ and the broadening matrix \mbox{\boldmath\Gamma}_{\alpha}. They are calculated as discussed in the previous section but the external bias is replaced by the KS bias and the Hxc gate is the long-time limit of the Hxc potential in the central region. We remind the reader that in the LB+DFT scheme the Hxc gate results from a self-consistent calculation of the density in the central region and there is no a priori reason for this Hxc gate to be the same as the long-time limit of the TDDFT Hxc gate.

It is worth emphasizing that the xc gate and bias in Eq. (76) are functionals of the density everywhere and at all previous times.

III.1 Linear response

In this section we work out and discuss the formula for the conductance resulting from the TDDFT formulation. To first order in the bias the current of Eq. (76) reads

[TABLE]

where

[TABLE]

is the transmission function calculated with equilibrium (zero bias) DFT Green’s function ${\mathbf{G}}$ and broadening matrix \mbox{\boldmath\Gamma}_{\alpha}. The conductance

[TABLE]

coincides with the KS conductance of the LB+DFT approach, see Eq. (21), only provided that $V_{\alpha,\rm xc}=0$. In the previous section we showed that the exact KS conductance and the exact conductance of the Anderson model are different in the Coulomb blockade regime. We therefore have to conclude that the xc bias is non zero in this case.

From linear-response TDDFT ewk.2004 ; skgr.2007

[TABLE]

where $f_{\rm xc}$ is the TDDFT kernel and $\delta n({\bf r},t)$ is the density variation. The assumption of a steady state implies that the kernel $f_{\rm xc}\rightarrow 0$ for $|t-t^{\prime}|\rightarrow\infty$ and that $\delta n({\bf r},t\rightarrow\infty)=\delta n_{\alpha}$ for ${\bf r}$ deep inside lead $\alpha$. In Eq. (81) the contribution of the molecular region to the spatial integral is negligible in the thermodynamic limit. Hence, it is convenient to define the quantity

[TABLE]

and rewrite Eq. (81) as

[TABLE]

where \mbox{\mathcal{S}}_{\beta}=\int_{\rm lead\;\beta}d{\bf r}^{\prime}_{\perp} is the area of the transverse section of lead $\beta$. Notice that due to charge conservation we have

[TABLE]

It could be tempting to gain some insight in the behavior of $f_{\rm xc}^{\alpha\beta}$ by performing equilibrium DFT calculations on leads of finite length and different densities. In doing so, however, we would get the equilibrium DFT kernel which corresponds to taking the limit $t\rightarrow\infty$ before the limit ${\mathbf{r}}\to\alpha$. In fact, the equilibrium DFT kernel is the response of the equilibrium xc potential to a density variation and, deep inside the leads, is determined by the condition of charge neutrality alone. As, in general, the limit $t\rightarrow\infty$ and ${\mathbf{r}}\to\alpha$ do not commute we cannot model $f_{\rm xc}^{\alpha\beta}$ using leads of finite length.

Inserting the expression for $\delta V^{\alpha}_{\rm xc}$ into Eq. (78) we find

[TABLE]

where

[TABLE]

The expression for $\delta I$ in Eq. (85) is correctly gauge invariant. Under a gauge transformation the kernel $f_{\rm xc}({\mathbf{r}},{\mathbf{r}}^{\prime})$ changes by the addition of an arbitrary function $q({\bf r})+q({\bf r}^{\prime})$ hg.2012 and $\Phi$ is invariant under this transformation. In conclusion

[TABLE]

The quantity v\equiv\delta I/(\mbox{\mathcal{S}}_{L}\delta n_{L}) is the speed of the charge wavefront propagating in the leads after the sudden switch-on of the external bias UimonenKhosraviStanStefanucciKurthLeeuwenGross:11 and it is of the order of the Fermi velocity. In the following we refer to $\Phi G_{s}/v$ as the dynamical xc correction since $\Phi$ is expressed in terms of the TDDFT kernel. An equation similar to Eq. (87) can also be obtained within the framework of time-dependent current density functional theory as has been shown in Ref. kbe.2006, .

III.2 Anderson Model in the Coulomb Blockade regime

We go back to the Anderson model at temperatures higher than the Kondo temperature but smaller than the charging energy $U$. This is the so called Coulomb blockade (CB) regime where the Abrikosov-Suhl resonance (or Kondo peak) in the spectral function has disappeared. However, as we already pointed out in Sec. II.4.1, the CB peaks present in the exact interacting conductance are completely absent in the KS conductance. We would like to emphasize that the physical situation discussed here is distinct from the one of Ref. tfsb.2005, . In Ref. tfsb.2005, the discontinuity is responsible for keeping the HOMO doubly occupied and the LUMO empty as the gate potential becomes more attractive (closed shell). In this case the discontinuity correctly suppresses $G_{s}$ at even $N$. Instead, at odd $N$ the discontinuity has the opposite effect since it pins the KS gate to the Fermi energy, thereby favouring the tunneling of electrons. Open-shell molecules in the CB regime represent a striking example of the inadequacy of standard DFT transport calculations. We now show that this is due to the lack of the dynamical xc correction discussed in Sec. III.1.

To gain some insight into the density dependence of $\Phi/v$ we reason as follows. Away from half-filling the MB and KS system behave similarly and consequently $G\simeq G_{s}$. On the other hand at half-filling, i.e., $N=1$, the exact conductance is strongly suppressed whereas the KS conductance is of the order of the quantum of conductance $G_{0}$. Therefore the dynamical xc correction has to be small for $N\neq 1$ and large for $N=1$. Interestingly, this is the same behavior of the derivative of the Hxc potential $\partial v_{\rm Hxc}/\partial N$. In fact, we can demonstrate that the two quantities are intimately related.

We consider the Anderson model in equilibrium and calculate the compressibility $\kappa=\partial N/\partial\mu$ using

[TABLE]

where $A$ is the interacting spectral function. It is a matter of simple algebra to show that

[TABLE]

where $G$ is the interacting conductance of Eq. (55) and in the last equality we defined

[TABLE]

The interacting and DFT compressibilities are the same by construction since the Hxc potential is such that the interacting and DFT densities are identical. Therefore we can also write

[TABLE]

where $G_{s}$ is the KS conductance and

[TABLE]

is the KS spectral function. The latter depends on $\mu$ through $N$ and the dependence on $N$ is all contained in $v_{\rm Hxc}$. Taking into account that $\frac{\partial A_{s}}{\partial v_{\rm Hxc}}=-\frac{\partial A_{s}}{\partial\omega}$ we have

[TABLE]

Inserting this result into Eq. (91), solving for $\kappa$ and equating the interacting and DFT expressions we find

[TABLE]

We observe that this relation is valid for any temperature; no approximations have been made so far.

We are interested in modelling the dependence of $R$ on $N$ for temperatures in the CB regime. In this case the many-body (MB) spectral function $A\simeq A^{\rm mod}$, see Eq. (58), and therefore $R(v)=I(v)-I(v+U)$ where $I(E)\equiv\int f(\omega)\ell_{\gamma}(\omega-E)$. The relation between $v$ and $N$ stems from Eq. (88) and reads

[TABLE]

from which it follows that

[TABLE]

For $v<\mu$, or equivalently for $N<1$, we have $I(v+U)\ll 1$. Thus for $N<1$ we can write $N\simeq 2I(v)/(1+I(v))$ and solving for $I(v)$ we get $I(v)\simeq N/(2-N)$. The expression of $I(v)$ for $N>1$ can be inferred using the ph symmetry and the final result is

[TABLE]

where $\delta N=N-1$. Inserting Eq. (97) into Eq. (94) we obtain the following DFT result for the conductance

[TABLE]

This relation is of great utility since it allows to estimate the dynamical xc correction from equilibrium DFT. In fact, the dynamical xc correction of Eq. (87) is entirely expressed in terms of equilibrium DFT quantities. Moreover, whereas $\Phi$ involves the TDDFT kernel with coordinates in the leads the correction in Eq. (98) involves only the DFT $v_{\rm Hxc}$ in the molecular junction. The accuracy of Eq. (98) is examined in Fig. 9, and benchmarked against the interacting conductance of Eq. (56). Even though the approximate $R$ is not on top of the exact one, see inset, the agreement between the two conductances is extremely good. Most importantly the plateau of $G_{s}$, see Fig. 4, is completely gone.

III.2.1 Seebeck coefficient of the Anderson model

The idea to construct the dynamical xc correction to the conductance is rather general and can be used to calculate other response quantities using equilibrium DFT. As has been shown in Ref. YPKSD.2016, , one such quantity is the Seebeck coefficient that, according to a recently proposed DFT framework for thermal transport EichVentraVignale:14 ; EichPrincipiVentraVignale:14 , does also contain dynamical xc corrections.

The Seebeck coefficient $S$ is defined as the ratio $S=(\delta V/\delta T)_{I=0}$, where $\delta V$ is the voltage that must be applied to cancel the current $\delta I$ generated by a small temperature difference $\delta T$ between the left and right leads. This definition corresponds to the phenomenological Seebeck coefficient of Refs. Mahan:06, ; Apertet:15, . For the Anderson model the Seebeck coefficient takes the form Dong:02

[TABLE]

with, according to our notation, $f^{\prime}\equiv df/d\omega$ and $A$ the interacting spectral function. To obtain an expression for $S$ in terms of equilibrium DFT quantities we calculate $dN/dT$ from Eq. (88). In the CB regime $A\simeq A^{\rm mod}$ depends on $T$ and $\mu$ exclusively through $N$. Using $dA/dT=(dA/dN)(dN/dT)$ it is easy to show that

[TABLE]

where $R$ is defined as in Eq. (90). Therefore, the numerator of the Seebeck coefficient in Eq. (99) is related to the temperature derivative of $N$. On the other hand, the denominator in Eq. (99) is related to the compressibility $\kappa$ of Eq. (89) since $G=-\frac{\gamma}{2}\int\!\frac{{\rm d}\omega}{2\pi}f^{\prime}(\omega)\,A(\omega)$, see Eq. (55). Therefore, we can write the Seebeck coefficient for the Anderson model in the CB regime as

[TABLE]

This is the expression we were looking for as both derivatives $dN/dT$ and $dN/d\mu$ can be calculated from equilibrium DFT. In the KS system $N=2\int\!\frac{{\rm d}\omega}{2\pi}f(\omega)A_{s}(\omega)$ where $A_{s}$ is the KS spectral function of Eq. (92). Since the Hxc potential $v_{\rm Hxc}$ depends on $N$ and $T$, at self-consistency $A_{s}$ depends implicitly (through $N$) on $\mu$ and both implicitly (through $N$) and explicitly on $T$. By calculating the required density derivatives and using $\frac{dv_{\rm Hxc}}{dT}=\left(\frac{\partial v_{\rm Hxc}}{\partial N}\right)_{T}\frac{dN}{dT}+\left(\frac{\partial v_{\rm Hxc}}{\partial T}\right)_{N}$, we obtain the relation

[TABLE]

In this result the KS Seebeck coefficient $S_{s}$ is defined as in Eq. (99) but with spectral funcion $A_{s}(\omega)$ in place of $A(\omega)$ [notice that the only requirement for the derivation of Eq. (102) is that $A_{s}$ is a function of $(\omega-v-v_{\rm Hxc})$]. The KS Seebeck coefficient is precisely the coefficient predicted by the LB+DFT approach which lacks the dynamical xc correction $\left(\frac{\partial v_{\rm Hxc}}{\partial T}\right)_{N}$. Although $v_{\rm Hxc}$ depends very weakly on temperature, it turns out that this weak dependence is still crucial to reproduce the MB Seebeck coefficient. Let us illustrate this point in more detail.

In the CB regime $\gamma$ is the smallest energy scale. We consider the limit of very weak contacts $\gamma\ll T,U$ and approximate $v_{\rm Hxc}$ by the exact Hxc potential of the isolated ($\gamma=0$) impurity sk.2011 , see Eq. (35). Having an analytic expression for $v_{\rm Hxc}$ we can evaluate both terms on the r.h.s. of Eq. (102). In Fig. 10 we show $S$ calculated from Eq. (99) (black) versus the gate $v$ and compare it with $S$ calculated from Eq. (99) using the MB spectral function of Eq. (58) (blue) as well as with $S$ calculated using the Rate Equation (RE) approach of Ref. BeenakkerStaring:92, (red), which is exact in the limit $\gamma\to 0$. All three approaches give the same Seebeck coefficient and densities (see inset).

It is instructive to analyze how the two terms in Eq. (102) contribute separately to yield to correct Seebeck coefficient. We see in Fig. 10 that the KS Seebeck coefficient $S_{s}$ (green) accounts for the correct linear behavior (with slope proportional to $T^{-1}$) at large values of $|v|$. This can easily be understood since as $\gamma\to 0$ the KS spectral function $A_{s}(\omega)\to 2\pi\delta(\omega-v-v_{\rm Hxc})$ and consequently $S_{s}=-(v+v_{\rm Hxc})/T$. Less obvious, instead, is the plateau of $S_{s}$ for $v\in(-U,0)$. The plateau is a direct consequence of the step in $v_{\rm Hxc}[N]$ at $N=1$, responsible for blocking electrons with energy below $v+U$ from entering the impurity site (see inset). The CB opens a gap in the noninteracting straight line $-v/T$, shifting it leftward by $U$ for $v<-U$ and generating the correct behavior at large negative values of $v$. We may say that the plateau is a manifestation of the CB, an equilibrium property of the impurity occupation. Although the KS Seebeck coefficient correctly captures the aforementioned gap it entirely misses the oscillation of $S$ for $N\approx 1$. The dramatic consequence of this fact is that $S_{s}$ severely underestimates the interacting Seebeck coefficient. It is remarkable that this problem is perfectly cured by the dynamical xc correction $\partial v_{\rm Hxc}/\partial T$, see Fig. 10 (cyan). Thus, the explicit temperature dependence of $v_{\rm Hxc}$ is the key ingredient for the Seebeck coefficient not to vanish in the CB regime BeenakkerStaring:92 ; StaringMolenkamp:93 ; Dzurak:93 ; Dzurak:97 .

One last remark before closing this section. The use of any temperature-dependent LDA potential GuptaRajagopal:82 ; ksdt.2014 in a LB+DFT calculation would not only miss the oscillation induced by $\partial v_{\rm Hxc}/\partial T$ but also the plateau in $S_{s}$ due to the lack of the step in $v_{\rm Hxc}$ at $N\approx 1$.

III.3 Constant Interaction Model in the Coulomb Blockade regime

In this section we extend the analysis on the Anderson model to the Constant Interaction Model (CIM) introduced in Sec. II.3.2. For simplicity we assume that the broadening matrix \mbox{\boldmath\Gamma}_{\alpha,mn}=(\gamma/2)\delta_{mn} is diagonal and proportional to the identiy matrix. In Sec. II.3.2 we proved that the Hxc potential of the zero-temperature isolated CIM is a uniform shift depending on the total number of particles only, i. e., $v_{{\rm Hxc}}[n]({\mathbf{r}})=v_{\rm Hxc}[N]$, see Eq. (43). For a broadening $\gamma$ and temperature $T$ much smaller than both the level spacings and the charging energies, the inhomogeneity of the Hxc potential as well as the dependence of the Hxc potential on the local occupations can be safely discarded. Then we can go through the same steps of the single-level derivation of Sec. III.2 and find again Eq. (98). The only difference is that $\delta N$ is given by the deviation of $(N-{\rm Int}[N])$ from unity.

To illustrate the importance of the dynamical xc correction to the conductance we approximate $v_{\rm Hxc}$ as in Eq. (67). We recall that the charging energies $U_{J}$ are given by the xc part of the derivative discontinuity of the CIM with $J$ electrons PerdewParrLevyBalduz:82 . For the widths we take $W_{J}=0.16\,\gamma/U_{J}$ which is consistent with Ref. ES.2011, . We mention that in the CB regime the temperature $T>\gamma$ and hence the smeared steps of $v_{\rm Hxc}$ should be broadened by $T$ and not by $\gamma$ as in Eq. (67). Nevertheless, as we already discussed, this quantitative feature has no effect in the qualitative behavior of the number of particles $N$ and KS conductance $G_{s}$ as functions of the gate $v$.

Possible physical realizations of the CIM are quantum dots made from metallic single-wall nanotubes (SWNT) lbp.2002 ; BBNIS.2002 ; SJKDKvdZ.2005 . This has been shown by Oreg et al. OBH.2000 who were able to reproduce the observed fourfold periodicity in the electron addition energy of a SWNT of finite length using the following CIM Hamiltonian

[TABLE]

Here $\sigma$ is the spin index, $\nu=0,1$ is the band index and $l$ is the integer of the quantized quasi-momentum of the electrons. The entire Hamiltonian is expressed solely in terms of the occupation numbers $n_{l\nu\sigma}$ since $N_{\sigma}\equiv\sum_{l\nu}n_{l\nu\sigma}$ (total number of electrons with spin $\sigma$) and $N=N_{\uparrow}+N_{\downarrow}$ (total number of electrons). The finite length of the SWNT causes a finite subband mismatch $\delta$ so that the single-particle energies are

[TABLE]

with $\Delta$ the average level spacing. In Eq. (103) $E_{C}$ is the charging energy ($N_{0}$ is the number of electrons of the charge neutral SWNT quantum dot), $\delta U$ is the extra charging energy for two electrons in the same energy level and $J$ is the exchange energy between electrons of opposite spin. With the parameters of Ref. OBH.2000, , when an extra electron enters the nanotube it occupies the lowest available single-particle energy level. Thus the aufbau is the same as that of the noninteracting Hamiltonian. Accordingly the spin of the ground state is [math], $1/2$, [math], $1/2$, …for 0, 1, 2, 3, …extra electrons. Let $E(N)$ be the ground state energy of the SWNT with $N$ extra electrons. Using Eq. (103) it is straightforward to obtain

[TABLE]

and so on. The function $E(N)$ with $N$ a real continuous variable has a discontinuous derivative at integers $N$ and the size of this discontinuity is given by $\Delta(N)=E(N+1)-2E(N)+E(N-1)$, see Ref. PerdewParrLevyBalduz:82, . One finds

[TABLE]

and $\Delta(N)=\Delta(J)$ if $N-J=0$ mod2 with $J=1,2$.

The Hxc potential $v_{\rm Hxc}[N]$ has the property that for a given chemical potential the ground state occupations $n_{l\nu\sigma}$ of the KS system are the same as those of the CIM. From Eqs. (106) we find that the charging energies $U_{J}$ are given by

[TABLE]

and $U_{J}=U_{K}$ if $J-K=0$ mod2 with $K=1,2$. The average values of these parameters can be found in Ref. OBH.2000, . In order to match the position of the conductance peaks of the SWNT of length $\simeq 100$ nm we use (all energies are in meV): $\Delta=9.2$, $\delta=2.27$, charging energy $E_{C}=2.485$, exchange energy $J_{\rm x}=0.7$, extra charging energy for doubly occupied levels $\delta U=0.37$ KS.2013 .

The total number of particles is obtained from the self-consistent solution of the KS equation (19) which in equilibrium can be written as

[TABLE]

where ${\mathbf{A}}_{s}(\omega)={\mathbf{A}}_{L,s}(\omega)+{\mathbf{A}}_{R,s}(\omega)$ is the total KS spectral function and ${\mathbf{G}}$ is the Green’s function of Eq. (69). Taking into account that the broadening matrices are diagonal and that $v_{\rm Hxc}$ is uniform then ${\mathbf{A}}_{s}$ is diagonal. For the SWNT with Hamiltonian in Eq. (103) the trace of the KS spectral function reads

[TABLE]

In Fig. 11 we plot $v_{\rm Hxc}[N]$ as well as $N$ at self-consistency as a function of $v_{g}=v_{0}+\alpha v$ for two different values of $\gamma$. The potential energy $v_{0}$ is determined by requiring that our reference energy is the same as in Ref. lbp.2002, . The parameter $\alpha=C/C_{g}$ is the ratio between the total capacitance and the gate capacitance. For the experiment in Ref. lbp.2002, this ratio is about 250.

To obtain the KS conductance at a certain value of $v$ we evaluate the KS spectral function with $N=N[v]$ and then calculate the KS conductance in accordance with Eq. (21), i.e.,

[TABLE]

Subsequently we correct $G_{s}$ according to

[TABLE]

see Eq. (98). The KS conductance $G_{s}$ as well as the conductance with dynamical xc corrections $G$ are shown in Fig. 12 for different values of the broadening parameter $\gamma$. As expected, for small $\gamma$ the KS conductance behaves like in Fig. 7, i. e., it exhibits a Kondo plateau whenever the number of particles $N$ is close to an odd integer. The dynamical xc correction suppresses this plateau and yields the correct CB pattern.

In Fig. 13 we compare $G_{s}$ and $G$ with the experimental conductance. The conductance $G$ represents a considerable improvement over $G_{s}$ which, instead, shows two deformed Kondo plateaus per period. Notice that the fourfold periodicity is also captured. More details can be found in Ref. KS.2013, .

III.3.1 Seebeck coefficient of the CIM

For temperatures $T\gg\gamma$ the Seebeck coefficient of the CIM exhibits a sawtooth behavior as a function of gate voltage, with “jumps” occurring when the number of particles crosses an integer. In addition to these jumps, a superimposed fine structure of wiggles spaced by $\Delta\varepsilon$ emerges whenever the level spacing $\Delta\varepsilon$ is larger than the temperature BeenakkerStaring:92 . The wiggles are associated to charged excitations from the ground state with $(N-1)$ particles to some excited state with $N$ particles.

For diagonal broadening matrices \mbox{\boldmath\Gamma}_{\alpha,mn}=\delta_{mn}\gamma/2 we can again express the Seebeck coefficient as in Eq. (101). The derivation is identical provided that we replace the interacting spectral function $A$ with its trace ${\rm Tr}[{\mathbf{A}}]$. Since $N$ can be calculated from DFT, Eq. (101) allows us to express the Seebeck coefficient in a pure DFT fashion. Approximating the Hxc potential as a uniform shift, see discussion at the beginning of Sec. III.3, it is straightforward to show that

[TABLE]

Like in Eq. (67) we construct $v_{\rm Hxc}[N]$ by summing over all possible charged states the Hxc potential of the Anderson model. However, due to the importance of the temperature dependence we use the finite-temperature Hxc potential of the isolated Anderson impurity, see Eq. (35). Thus we have

[TABLE]

where $U_{J}$ is the charging energy and the extended $g^{\rm ext}_{U}$ function is defined according to

[TABLE]

with $g_{U}$ given in Eq. (36). Like the Hxc potential in Eq. (67) also this Hxc potential has a staircase behavior with smeared steps of width $U_{J}$ between two consecutive integers but the smearing is governed by $T$ instead of $\gamma$.

It is worth noting that the property of the Hxc potential of being the same for all energy levels is an exact feature only at zero temperature. In fact, our approximate $v_{\rm Hxc}$ can reproduce only the occupations corresponding to a thermal mixture of ground states with different number of particles. To illustrate what physics is lost in this way we first consider a two-level CIM with repulsion energy $U$. We choose a temperature $T=0.03$ much larger than $\gamma=0.001$ and the energy of the levels $\varepsilon_{i}=\varepsilon_{i}^{0}+v$ with $\varepsilon^{0}_{1}=0$ and $\varepsilon_{2}^{0}=0.3$ and $v$ the external gate potential. Here all energies are given in units of $U$. The left panel of Fig. 14 shows the total number of particles $N=n_{1}+n_{2}$ as well as the occupation $n_{2}=\sum_{\sigma}n_{2\sigma}$ of the highest level as obtained using DFT with Hxc potential in Eq. (113) and the RE approach. The approximation of a uniform Hxc potential has no effect on $N$, which is identical in both approaches, but it introduces exponentially small discrepancies in $n_{2}$ (and hence in $n_{1}$). These discrepancies are due to the neglect of excited states in the thermal mixture represented by $v_{\rm Hxc}$. Therefore, we expect that some of the wiggles in the Seebeck coefficient are not captured by our approximation (which accounts only for the addition of electrons in the lowest available level). This is confirmed by the right panel of Fig. 14 where the wiggles associated to the addition energies of excited states emerge using a rate equations approach (red) and are absent using our approximate DFT treatment (black). Nevertheless, we emphasize that the wiggles stem from $S_{s}$ and are not due to the xc correction. The latter is responsible for the large sawtooth oscillations and, as Fig. 14 clearly shows, it is the dominant contribution to $S$.

Experimental measurements of the Seebeck coefficient for an individual single-wall carbon nanotube in the CB regime have been reported in Ref. SPK.2003, . In order to show the performance of our DFT scheme we extracted both single-particle energies and charging energies from the experimental results. We again consider the Hxc potential of Eq. (113) but, in contrast to the model calculations described previously, the charging energies $U_{J}$ depend on the charging state $J$. Details on the parameters can be found in Ref. YPKSD.2016,

In Fig. 15 we show the interacting conductance $G$ calculated using Eq. (111), see also Ref. KS.2013, , (upper panel) and the Seebeck coefficient $S$ calculated from Eq. (112) (lower panel) versus the gate voltage $v$ for temperature $T=4.5$ K and coupling $\gamma=0.02$ meV. For comparison we also show the KS Seebeck coefficient $S_{s}$. The latter fails in reproducing the characteristic sawtooth behaviour of the experimental results. Instead, the interacting Seebeck coefficient calculated from Eq. (112) clearly shows the peak and valley structures observed in experiment, confirming again the crucial role of the xc correction. Remarkably, the fine structure wiggles (kinks in some cases) are correctly captured too.

IV Steady-State Density Functional Theory for Transport at Finite Bias

In Sec. II.1 we have discussed the standard approach to transport within a DFT framework, the LB+DFT formalism which combines static (ground state or equilibrium) DFT with the Landauer-Büttiker formalism. We have pointed out that formally this approach is incomplete as even the knowledge of the exact xc potential of DFT doesn’t guarantee that the corresponding (steady-state) currents are exact. In the regime of linear response we have shown with explicit examples that LB+DFT transport properties such as the zero-bias conductance or the Seebeck coefficient can capture some but not all of the correct physics of transport.

In Sec. III, on the other hand, we have discussed TDDFT as a truly non-equilibrium density functional approach which in principle can describe electronic transport correctly. This is of course true in the time domain if one is interested in explicitly time-dependent currents such as transients or currents in response to a time-dependent bias. However, it is also true for the steady-state regime which is interpreted as the long-time limit of the time evolution of a system driven out of equilibrium by a DC bias. In fact, we have seen in Eq. (76) that TDDFT leads to an expression for the steady current which is structurally identical to the original LB+DFT expression (20) but with the crucial difference that the applied bias has to be corrected by an xc contribution.

In principle, the TDDFT xc potential (and thus also the xc correction to the bias) is a functional with “memory”, i.e., it depends not only on the instantaneous density but also on its entire history at previous times. Moreover, it not only depends on the density in the central device region but also on the one deep inside the leads. On the other hand, for the models we studied we have derived expressions for the xc correction to the bias (at least in the linear regime) which depend only on the density of the device. If we are honest, though, while our interpretation of the xc bias correction was based on TDDFT ideas, its derivation was not. This latter fact maybe shouldn’t come as a surprise: the construction of TDDFT functionals beyond the adiabatic approximation is a notoriously difficult task.

In the present Section we will present a DFT formulation, which we will call i-DFT, for steady-state transport, i.e., we here aim to describe a system in its (non-equilibrium) steady state and are not interested in the time evolution towards this state (just as the LB+DFT formalism does). We will try to incorporate some of the lessons learned from previous (model) studies into the new framework, one of which is the overriding importance of an xc correction to the bias, another one the formulation in terms of quantities defined in the device region only. At the very fundamental level we will first choose a set of basic “density” variables for which, under certain conditions, one can prove a one-to-one correspondence between this set of densities and a set of potentials. Once the basic formalism is established, we will again focus on its usefulness in the description of transport through strongly correlated systems. First we will show how the new formalism can handle Coulomb blockade both at zero and at finite bias. Then we will apply the insights gained here again to the problem of the SIAM and will construct a functional which is able to correctly describe both the Kondo and the Coulomb blockade regime as well as the transition from one to the other.

IV.1 Foundations of the i-DFT formalism

As before, we consider a central region attached to left and right leads. In equilibrium, the system is subject to an electrostatic potential $v_{0}({\bf r})$, e.g., the potential generated by the nuclei. Out of equilibrium, the system is exposed to an external bias potential $v_{b}({\bf r})$ generated by an external battery. In addition the system may also be subject to an additional gate potential $v_{g}({\bf r})$ which vanishes deep inside the leads. In these regions of space the classical potential $v_{0}({\bf r})+v_{b}({\bf r})+v_{\rm H}({\bf r})$ (where $v_{\rm H}({\bf r})$ is the Hartree potential), differs by a uniform shift $V/2$ ($-V/2$) from its equilibrium value, where $V$ is the potential drop (bias) across the junction. Now we split the different components of the potentials according to the different regions, e.g., we write $v_{0}({\bf r})=v_{0,L}({\bf r})+v_{0,C}({\bf r})+v_{0,R}({\bf r})$ where $v_{0,\alpha}({\bf r})=v_{0}({\bf r})$ for ${\bf r}\in\alpha$ and zero otherwise. The total potential in region $C$ then is $v_{C}({\bf r})=v_{0,C}({\bf r})+v_{{b},C}({\bf r})$. We can also split the density in a similar way, i.e., $n({\bf r})=n_{L}({\bf r})+n_{C}({\bf r})+n_{R}({\bf r})$. With these definitions we will show below that there is a one-to-one correspondence between the pair $(v_{C}({\bf r}),V)$ and the pair $(n_{C}({\bf r}),I)$ where $I$ is the steady-state current through region $C$. The fundamental theorem of i-DFT can then be formulated as

Theorem: For any finite temperature and for fixed potentials $v_{\alpha}({\bf r})$ in leads $\alpha\in\{L,R\}$ the map $(v_{C}({\bf r}),V)\longrightarrow(n_{C}({\bf r}),I)$ is invertible in a finite bias window around $V=0$.

Proof: In order to prove the theorem we show that the Jacobian

[TABLE]

is non-vanishing. The upper left block $\chi_{C}({\bf r},{\bf r}^{\prime})=\frac{\delta n_{C}({\bf r})}{\delta v_{C}({\bf r}^{\prime})}_{V=0}$ is the static, equilibrium density response function for the contacted $L-C-R$ system but evaluated for both ${\bf r}$ and ${\bf r}^{\prime}$ in region $C$ and $G\equiv\frac{\partial I}{\partial V}$ is the zero-bias conductance. As a first observation we note that the variation \frac{\delta I}{\delta v_{C}({\bf r}^{\prime})}\big{|}_{V=0} vanishes since at zero bias a change in the central potential does not induce a steady current. Thus, it remains to be shown that both entries on the diagonal of the Jacobian (115) have a definite sign. We first look at the equilibrium density response function $\chi_{C}({\bf r},{\bf r}^{\prime})$ which can be calculated using leads of finite length $L$ at taking the limit $L\to\infty$ at the end. At finite temperature $1/\beta$ and at chemical potential $\mu$, the Lehmann representation of $\chi_{C}({\bf r},{\bf r}^{\prime})$ reads

[TABLE]

Here $Z$ is the partition function and the sum is over a complete set of many-body eigenstates $|\Psi_{i}\rangle$ of the contacted system with energy $E_{i}$ and particle number $N_{i}$. Furthermore we have defined the excitation energies $\Omega_{ij}=E_{i}-E_{j}$, the excitation amplitudes $f_{ij}({\bf r})=\langle\Psi_{i}|\hat{n}({\bf r})|\Psi_{j}\rangle-\delta_{ij}n({\bf r})$ (with the density operator $\hat{n}({\bf r})$) and $\eta$ is a positive infinitesimal. In order to prove the invertibility of $\chi_{C}({\bf r},{\bf r}^{\prime})$ we have to show that for an arbitrary test function $t({\bf r})$ we have

[TABLE]

with $T_{ij}\equiv\int_{C}{\rm d}^{3}r\;f_{ij}({\bf r})t({\bf r})$. It is easy to see that for $E_{i}\neq E_{j}$ we have $\Omega_{ij}\left(e^{-\beta E_{i}}-e^{-\beta E_{j}}\right)<0$ and thus the l.h.s. of Eq. (117) can be zero only if $T_{ij}=0$ for any $i,j$ with $E_{i}\neq E_{j}$. Obviously, for an arbitrary test function $t({\bf r})$ this cannot happen StefanucciKurth:15 and we have to conclude that $\chi_{C}$ is invertible.

Also for the zero-bias conductance one can write down the Lehmann representation bsw.2006 which reads

[TABLE]

where $I_{ij}\equiv\langle\Psi_{i}|\hat{I}|\Psi_{j}\rangle$ with the longitudinal current operator $\hat{I}$. From this expression one can see that $G>0$. Therefore we find that the Jacobian $J_{V=0}=\det[\chi_{C}]G<0$. Since $J_{V}$ is a continuous function of $V$ around $V=0$, there exists a finite interval around $V=0$ (whose size depends on $v_{C}$) for which $J_{V}<0$. Therefore in this domain the map $(v_{C}({\bf r}),V)\longrightarrow(n_{C}({\bf r}),I)$ is invertible.

In what follows, we will omit the subscript $C$ again but it is understood that all local quantities (potentials, densities) refer to the central device region only. Let $(n({\bf r}),I)$ be the density and steady current induced by the potentials $(v({\bf r}),V)$ in an interacting junction. We assume that the pair $(n({\bf r}),I)$ is non-interacting $v$-representable, i.e., it is a physically realizable pair for a non-interacting system as well. Then our theorem guarantees that the pair of potentials $(v_{s}({\bf r}),V_{s})$, which leads to the same density and current in a non-interacting system, is unique. Following the usual KS procedure we can then define the Hxc gate potential and the xc bias as

[TABLE]

[TABLE]

Of course, these are purely formal definitions and in practice $v_{\rm Hxc}[n,I]$ and $V_{\rm xc}[n,I]$ have to be approximated.

The self-consistent KS equations of the i-DFT formalism are then given by Eqs. (16) and (20) with $V_{\alpha}\to V_{\alpha}+V_{\alpha,\rm xc}$. For a symmetric bias these equations read

[TABLE]

[TABLE]

where $s_{R/L}=\pm 1$. Although the i-DFT equations are very similar to the LB+DFT equations there are important differences: in LB+DFT, the Hxc gate potential is a functional of the density alone and therefore there is a self-consistency condition only for the density, see Eq. (16). The current is then evaluated with this self-consistent Hxc gate potential from Eq. (20) which has the same structure as Eq. (122) but with vanishing xc bias. In i-DFT, on the other hand, we have to take into account the generally non-vanishing xc contribution to the bias. Both this xc bias as well as the Hxc gate depend on density and current and thus the self-consistency conditions for these two quantities are coupled and have to be solved together.

Also in comparison to TDDFT (for an applied DC bias in the long-time limit), the i-DFT equation for the current is structurally identical to Eq. (76) sa-1.2004 ; sa-2.2004 . However, the TDDFT Hxc gate potential and xc bias are functionals of the density everywhere, i.e., both in the device region and in the leads. Furthermore, the TDDFT potentials at given time $t$ depends on the full history of the density at all previous times. In contrast, the i-DFT functionals are independent of history (which is not surprising since i-DFT only deals with steady states) and depend only on the density in the device region as well as on the steady current. The augmented local character of the i-DFT xc potentials agrees with similar findings in time-dependent current density functional theory v.1995 ; nptvc.2007 .

The zero-bias conductance in the i-DFT formalism can be derived by linearising Eq. (122) in the bias leading to the simple but exact result StefanucciKurth:15

[TABLE]

Compared to the TDDFT result (87) for the conductance, which involves the zero-frequency and zero-momentum limit of the xc kernel, the above expression is more transparent. Its simple form will later also give a hint on the design of approximate i-DFT functionals.

IV.2 i-DFT functionals for the Coulomb blockade regime

So far we have presented the formal foundations of the i-DFT framework. In order for i-DFT to be applied, however, we need approximate xc functionals. In the present Section we will construct such approximations for the model systems we studied before, i.e., the SIAM and the CIM. For now, our aim is to construct approximations which capture the essential physics of Coulomb blockade both at zero and at finite bias. The construction of these approximations will be done by reverse-engineering from standard techniques typically used to describe Coulomb blockade.

We start again with the SIAM and again, for simplicity, we restrict ourselves to the wide-band limit. In order to obtain the i-DFT xc potentials we need a model for the density and current of the biased, interacting system. Fortunately, such a model is easily constructed using the ingredients already introduced in Sec. II.4.1. Using the model spectral function $A^{\rm mod}(\omega)$ of Eq. (58), the steady-state density $N=n$ and current $I$ of the SIAM can be calculated from

[TABLE]

and

[TABLE]

The resulting densities and currents are in excellent agreement with the results of the rate equations rate-paper1 ; rate-paper2 , the standard technique to describe Coulomb blockade for weakly coupled systems. In the reverse-engineering procedure we numerically invert Eqs. (124) and (125) for a given, fixed pair $(N,I)$ both for the interacting and the non-interacting system and then extract the i-DFT xc potentials according to Eqs. (119) and (120). One can actually prove StefanucciKurth:15 that for the SIAM in the WBL the map $(v,V)\to(N,I)$ is invertible for any value of the bias $V$ (infinite bias window). The domain spanned by $N$ and $I$ is $|I|\leq(\gamma/2)N$ for $N\in[0,1]$ and $|I|\leq(\gamma/2)(2-N)$ for $N\in[1,2]$.

The Hxc gate $v_{\rm Hxc}[N,I]$ and the xc bias $V_{\rm xc}[N,I]$ resulting from the reverse engineering procedure are shown in Fig. 16. The most prominent features are smeared steps of height $U/2$ for $v_{\rm Hxc}$ and of height $U$ for $V_{\rm xc}$ along the lines $N=1\mp I/\gamma$. The DFT xc discontinuity of $v_{\rm Hxc}[N,0]$ bifurcates as current starts flowing. The sign of the xc bias is opposite to the current, i.e., the effective KS bias $V+V_{\rm xc}$ is lower than the external bias $V$. This is in agreement with the model study of Ref. sds.2013, . The derivative $(\partial V_{\rm xc}/\partial I)_{I=0}<0$ and therefore, according to Eq. (123), $G_{s}$ becomes the upper limit of the interacting zero-bias conductance.

The i-DFT xc potentials at finite current can be parametrized in the spirit of Eq. (61) for the zero-current case as

[TABLE]

and

[TABLE]

where $W$ is defined according to Eq. (62) and we have introduced for later use an extra parameter which here we set to unity, $\lambda_{1}=1$. Note that Eq. (IV.2) reduces to Eq. (61) in the limit of zero current while Eq. (127) vanishes in this limit. We have verified that the self-consistent i-DFT results using the xc potentials of Eq. (IV.2) and (127) are in excellent agreement with the results of the rate equations (as they should be).

The reverse-engineering procedure can also be applied to the CIM. We consider a CIM with $M$ levels described by the Hamiltonian (40) and coupled to wide band leads, $\Gamma_{\alpha,ij}(\omega)=\delta_{ij}\gamma/2$. In general, the i-DFT xc potentials depend on all level occupations. For simplicity, we here restrict ourselves to the case of $M$ degenerate single-particle levels. In this case, by symmetry, the xc potentials become functionals only of the total occupation $N=\sum_{i=1}^{M}\sum_{\sigma}n_{i\sigma}$ and the total current $I$. Above the Kondo temperature $T_{\rm K}$, both $N$ and $I$ can be obtained by solving the rate equations rate-paper1 . For given $(N,I)$ we numerically invert the map $(v,V)\to(N,I)$ both for the interacting and the non-interacting case and then obtain the xc potentials according to Eqs. (119) and (120). Again, in the wide-band limit the map is invertible for all $V$ and the codomain is $|I|\leq(\gamma/2)N$ for $N\in[0,M]$ and $|I|\leq(\gamma/2)(2M-N)$ for $N\in[M,2M]$. In Fig. 17 we show the xc potentials for $M=3$. As in the SIAM, the Hxc gate (xc bias) potential exhibits smeared steps of height $U/2$ ($U$) with a rather complex pattern for the edges which follow piecewise straight lines in the $(N,I)$-plane. Again, the xc discontinuity at integer $N$ and $I=0$ bifurcates as the current starts flowing with the edges having different slopes depending on $N$. The edges connect “vertices”, i.e., points in the $(N,I)$-plane where two edges meet, and typically the edges change slope at the vertices. We denote by $\Delta_{K}^{(s)}(N,I)$ the piecewise linear function of $N$ and $I$ which vanishes along the step edge passing through $(K,0)$ with positive ($s=+1$) or negative ($s=-1$) slopes (see top panel of Fig. 17 for examples).

In order to model the reverse-engineered xc potentials we have to understand the position of the vertices. We realized an interesting duality StefanucciKurth:15 : the vertices occur exactly at those points in the $(N,I)$-plane which correspond to plateau values of the particle number and the current in the gate-bias $(v,V)$ plane. Moreover, these plateau values for $N$ and $I$ can be calculated for any degenerate $M$-level CIM using simple expressions. From the rate equations, the degenerate CIM with $M$ levels leads to $(2M+1)^{2}$ distinct plateau values. A plateau can be uniquely identified by a pair of two integers $(m,n)$ with $m,n=0,\ldots,2M$. In the $(m,n)$ plateau for $n\leq m$, the probability $P(q)$, $q=m,\ldots,n$ of finding $q$ particles are all identical and given by P^{-1}(q)\equiv P_{n\leq m}^{-1}=\sum_{j=m}^{n}\left(\!\!\begin{array}[]{c}2\mathcal{M}\\ j\end{array}\!\!\right), whereas all other probabilities vanish. The density and current of the $(m,n)$-plateau are then given by

[TABLE]

and

[TABLE]

Knowing the plateau values one can define the step edges $\Delta_{K}^{(s)}(N,I)$ and with these one can parametrize the Hxc gate and xc bias of the degenerate $M$-level CIM as

[TABLE]

and

[TABLE]

where, again, $W$ is defined according to Eq. (62) and $\lambda_{1}=1$. Again, for an $M$-fold degenerate CIM the self-consistent i-DFT results using the xc potentials (130) and (131) are in excellent agreement with the rate equation results.

So far, we have used the rate equations to construct i-DFT xc potentials for the degenerate case. For the non-degenerate CIM, according to our i-DFT philosophy, the xc potentials are now functionals of the local occupations $n_{i}$ and the current $I$ instead of the total $N$ and $I$. However, from what we have learned so far we can still construct useful approximations to the xc potentials for the non-degenerate case without having to do the full reverse engineering from the rate equations.

Let $n=\{n_{1},\ldots,n_{M}\}$ be the occupations of the levels ${1,\ldots,M}$ of an $M$-level CIM with arbitrary single-particle level structure. Let $\mathcal{M}_{p}[n]$ be the degeneracy of the $p$-th largest occupation and $\mathcal{D}[n]$ the number of distinct densities. For instance if $M=5$ and $n=\{\frac{1}{3},\frac{1}{2},\frac{1}{2},\frac{1}{3},\frac{1}{3}\}$ then $\mathcal{M}_{1}=2$, $\mathcal{M}_{2}=3$ and $\mathcal{D}=2$. We further define $\mathcal{N}_{p}[n]=2\sum_{q=1}^{p-1}\mathcal{M}_{q}[n]$ as the maximum number of particles in the first $(p-1)$ levels with degenerate occupations ($\mathcal{N}_{1}=0$). The degeneracies $\mathcal{M}_{p}$ are used to construct the following i-DFT potentials

[TABLE]

[TABLE]

The dependence on the local occupations enters exclusively through the $\mathcal{M}_{p}$. At the joining points ($N=\mathcal{N}_{p+1}[n]$ and $I=0$) between two consecutive $v_{\rm Hxc}^{(\mathcal{M}_{p}[n])}$ we add a discontinuity with slopes $\pm 2/\gamma$. In fact, the slope of the lines delimiting the domain of the i-DFT potentials of a $M$-fold degenerate CIM are independent of $M$ (see Figs. 16 and 17).

In order to show the performance of i-DFT we have calculated the finite-bias differential conductances of a benzene junction. Here we model the benzene molecule by a six-level CIM with $U=0.5$ eV and the single-particle energies $\varepsilon_{i}=\varepsilon_{i}^{0}+v$ with $\varepsilon_{1}^{0}=-\varepsilon_{6}^{0}=5.08$ eV, $\varepsilon_{2}^{0}=\varepsilon_{3}^{0}=-\varepsilon_{4}^{0}=-\varepsilon_{5}^{0}=2.54$ eV. These parameters are taken from a Pariser-Parr-Pople model of benzene BursillCastletonBarford:98 while the coupling to the leads is $\gamma=0.05$ eV. We show the differential conductances from LB+DFT (top panel), i-DFT (middle), and the rate equations (bottom). Due to the step structure in the (current-independent) Hxc gate potential, in LB+DFT we see Kondo plateaus at zero bias for odd occupations. Moreover, we also see “Kondo-like” plateaus at finite bias, a feature which is certainly unphysical. Furthermore, many of the lines with finite $dI/dV$ at finite bias which are present in the rate equations, in LB+DFT are simply missing. In i-DFT, on the other hand, we don’t see any Kondo plateaus neither at zero nor at finite bias. At zero bias, this is not surprising since we reverse-engineered our (H)xc potentials from the rate equations which doesn’t describe the Kondo effect. At finite bias it is comforting that i-DFT cures the spurious Kondo plateaus of LB+DFT. Furthermore, in the low-bias regime (delimited by the dashed horizontal lines) i-DFT describes all lines in the $dI/dV$ correctly as in the rate equations. At higher bias, on the other hand, some lines are missing. This, however, can be attributed to our particular approximation in Eq. (LABEL:xcgate_cim_gen) where the Hxc potential acts as a uniform constant shift for all levels.

IV.3 i-DFT functional for the SIAM: from Kondo to Coulomb blockade regime

In Section II.4 we have seen that, at zero temperature, already the LB+DFT formalism can correctly describe the Kondo plateau in the zero-bias conductance. In Section III, on the other hand, we have designed approximations to linear transport coefficients (zero-bias conductance, Seebeck coefficient) which capture the correct physics in the Coulomb blckade regime, assigning the corrections over LB+DFT to dynamical xc corrections of TDDFT. In the previous Section, in yet another DFT framework, we have constructed approximations to the i-DFT xc potentials which are designed to work in the Coulomb blockade regime, now not only at zero but also at finite bias. The question then arises if it is possible to design i-DFT functionals which capture both the Kondo and the Coulomb blockade regimes correctly and also the transition from one regime to the other. In the present Section we address this question for the SIAM KurthStefanucci:16 .

The KS self-consistency conditions in Eqs. (121) and (122) applied to the SIAM simplify to

[TABLE]

[TABLE]

where $s_{R/L}=\pm$ and the KS spectral function is given in Eq. (92).

At first we consider the case of zero temperature, $T=0$. We know from Sec. II.4.1 that if we use an approximate Hxc gate potential at zero current $v_{\rm Hxc}[N,I=0]$ which has a step at half occupation $N=1$, the resulting KS zero-bias conductance $G_{s}$ exhibits a Kondo plateau as function of gate voltage. This is exactly the case if we use the approximation of Eq. (IV.2) (which reduces to Eq. (61) for $I=0$). On the other hand, we also know that in order to obtain the interacting zero-bias conductance $G$, we have to correct $G_{s}$ according to Eq. (123). If we use Eq. (127) as an approximation to the xc bias, it is easy to show that $\partial V_{\rm xc}/\partial I|_{I=0}$ (and thus the correction to $G_{s}$) is non-vanishing. The resulting $G$ will, instead of showing the Kondo plateau, exhibit the typical two-peak structure associated with Coulomb blockade. It is thus immediately clear what we have to do in order to recover the Kondo plateau in the i-DFT framework: we have to design the approximation to the xc bias in such a way that at $T=0$ the correction to $G_{s}$ according to Eq. (123) vanishes. As a second requirement on our improved functional we want the Hxc gate at zero current to be not just qualitatively correct but also quantitatively as accurate as possible. Fortunately, in Ref. blbs.2012, an accurate, ready-to-use parametrization $v_{\rm Hxc}^{(0)}[n]$ of the Hxc gate potential at $T=0$ has been designed. We thus can make the following ansatz for our modified SIAM xc potentials KurthStefanucci:16

[TABLE]

and

[TABLE]

where $v_{\rm Hxc}[N,I]$ and $V_{\rm xc}[N,I]$ are the functionals of Eqs. (IV.2) and (127), respectively, with the parameter $\lambda_{1}$ introduced there now kept as a fitting parameter to be determined. For the choice of the function $\tilde{a}^{(0)}[I]$ there are a few restrictions: by symmetry, it should be an even function of the current and for the correction $\partial V_{\rm xc}^{\rm SIAM}\partial I|_{I=0}$ to vanish, its value at $I=0$ should be unity. Furthermore, we want the effect of $\tilde{a}^{(0)}$ to fade out for increasing currents since Eqs. (IV.2) and (127) already captured the high-current Coulomb blockade physics correctly. We thus suggest the following form:

[TABLE]

Finally, we fix the value of the parameter $\lambda_{1}=2$ such that the i-DFT $I-V$ characteristic at the particle-hole symmetric point is in good agreement with results from the functional renormalization group (fRG) EckelHeidrichJakobsThorwartPletyokhovEgger:10 .

In the left panel we compare $I-V$ characteristics from i-DFT with those of fRG for $v=-U/2$. The agreement is excellent with small deviations at small biases. These deviations become somewhat more pronounced as $U$ increases. In the right panel we compare again $I-V$ characteristics but this time at fixed value of $U/\gamma=1$ for different gates $v$. Again, the agreement is excellent. Note that the fitting parameter $\lambda_{1}$ was chosen only to optimize the agreement at $v=-U/2$ but we obtain excellent results also for other values of the gate.

After having fixed our parametrization for the SIAM (H)xc potentials at $T=0$, we now turn our attention to an extension for finite temperatures. Here we do not repeat all the details entering into the finite-$T$ approximations (which can be found in Ref. KurthStefanucci:16, ) but rather sketch the main physical ingredients entering its construction.

As we have already discussed in Sec. II.4.1, at the particle-hole symmetric point the zero-bias conductance $G^{\rm ph}$ of the SIAM at finite temperature is a universal function of $T/T_{\rm K}$ where $T_{\rm K}$ is the Kondo temperature defined in Eq. (54) Costi:00 ; AleinerBrouwerGlazman:02 . This universal function is shown in Fig. 6. When constructing our finite-$T$ approximations for the i-DFT functionals, we want to recover this exact property. This can be achieved by replacing $\tilde{a}^{(0)}[I]$ in Eqs. (136) and (137) by

[TABLE]

where $b^{(T)}[N]$ is chosen such that at $N=1$ the exact $G^{\rm ph}$ is recovered. The function $\tilde{a}^{(T)}[I]$ is obtained by just replacing the broadening $W$ in $\tilde{a}^{(0)}[I]$, see Eq. (138), by a temperature-dependent function $W(T)$ KurthStefanucci:16 defined as

[TABLE]

The same replacement $W\to W(T)$ is also done in the functional forms of $v_{\rm Hxc}[N,I]$ and $V_{\rm xc}[N,I]$. Physically, this replacement reflects the expectation that at small temperature the broadening in the steps of the (H)xc potentials is dominated by $\gamma$ while at large temperatures it is dominated by $T$. The particular form for $W(T)$ was chosen such as to best reproduce the fRG differential conductances of Ref. JakobsPletyukhovSchoeller:10, , see Fig. 20. The two physical ingredients for designing a finite-temperature approximations (reproduction of the universal $G^{\rm ph}$ as function of $T$, temperature dependent broadening of the step features) are augmented by some smaller tweaks to reproduce well KurthStefanucci:16 the fRG zero-bias conductances of Ref. JakobsPletyukhovSchoeller:10, .

In Fig. 20 a comparison of fRG and i-DFT finite-bias differential conductances at $v=-U/2$ is shown for two different interaction strengths and temperatures. By construction, the proper reduction of the Kondo peak at $V=0$ with increasing temperature is correctly reproduced. Also the Hubbard sidebands are reproduced with good accuracy.

Finally, in Fig. 21 we compare NRG and i-DFT zero-bias conductances as function of gate for two interaction strengths and various temperatures. In general, the agreement is very good. Only for high $U/\gamma$ (right panel) and low temperatures the shape of the side peaks is slightly different.

In summary we can say that we have indeed been able to design i-DFT functionals which capture both Kondo and Coulomb blockade physics in the SIAM as well as the transition between the two regimes as temperature increases. These functionals allow for the accurate calculation of densities and currents of the SIAM in the steady state over a wide range of parameters at negligible numerical cost.

Of course, the design of our functionals relied heavily on the availability of accurate solutions of the SIAM obtained from other methods, just like in standard DFT the construction of the local density approximation (LDA) relies on the availability of accurate xc energies of the uniform electron gas from Quantum Monte Carlo calculations. Nevertheless, with the design of our functionals we have explicitly demonstrated that it is indeed possible to accurately describe transport through strongly correlated systems with DFT, thus disproving common wisdom that DFT is not suited to deal with strong correlations.

At this point it is worth to emphasize that the construction of our i-DFT xc potentials was guided by a few rather simple considerations. The lessons learned here are easily transferrable to more complicated systems. For instance, the incorporation of Kondo physics in the i-DFT description of the contacted CIM can easily be achieved by modifying the CIM functionals of the previous Sec. IV.2 along similar lines to the ones discussed here for the SIAM: since the zero-bias KS conductance $G_{s}$ already contains the correct Kondo features also for multi-level systems (see Sec. II.4.2) one only has to ensure that the correction to $G_{s}$ vanishes for $I=0$. For instance, for a degenerate CIM with $M$ levels this can be achieved by using the following Hxc potentials (at $T=0$)

[TABLE]

[TABLE]

where $v_{\rm Hxc}^{(M)}$ and $V_{\rm Hxc}^{(M)}$ are given by Eqs. (130) and (131), respectively, with the only difference that we now use the parameter value $\lambda_{1}=2$. The function $\tilde{a}^{(0)}[I]$ is the same one as defined in Eq. (138) and for the equilibrium Hxc potential $\bar{v}_{\rm Hxc}^{(0)}$ we use

[TABLE]

Here we have defined the extended function

[TABLE]

with the parametrization $v_{\rm Hxc}^{(0)}[N]$ of the equilibrium SIAM Hxc potential of Ref. blbs.2012, . Note that for $M=1$, Eqs. (141) and (142) reduce exactly to Eqs. (136) and (137), i.e., our accurate $T=0$ parametrizations for the i-DFT (H)xc potentials of the SIAM.

In Fig. 22 we show the zero-temperature differential conductance in the gate-voltage plane of a degenerate HOMO-LUMO CIM for $U/\gamma=8$ (top panel) and $U/\gamma=4$ (bottom panel). One can clearly appreciate the Kondo strip at zero voltage for gates $v\in(-3U,0)$ as well as the fact that the height of the strip for $v\in(-2U,-U)$ is twice as large as the height for $v\in(-3U,-2U)$ and $v\in(-U,0)$. This simple incorporation of Kondo physics in more complicated models seems hard to achieve within other frameworks to deal with transport through correlated systems.

V Conclusion and outlook

The chemical complexity of molecular junctions calls for a first-principle description of the molecule attached to leads in order to make quantitative predictions and/or comparisons with experiments. DFT is a computationally efficient theory which, in principle, is well suited to achieve this goal. However, DFT strongly relies on the quality of the xc functional and, at present, most available approximations are inadequate for strongly correlated junctions. Furthermore, the applicability of (equilibrium) DFT in quantum transport needs to be carefully discussed.

In this Topical Review we have revisited the standard LB+DFT approach to quantum transport and showed that, even with the exact xc functional this approach is not exact. We have presented two exact frameworks which can cure the deficiencies of LB+DFT. The first is based on TDDFT and leads to an xc corrections to the applied bias. Unfortunately this xc bias correction as well as the Hxc potential in the molecular region are functionals of the density everywhere and at all previous times. This circumstance makes the development of practical approximations very difficult, especially if we are not interested in the full time evolution but only in steady-state properties. The second theory is i-DFT and, like TDDFT, predicts the existence of an xc bias correction. However, the basic variables of i-DFT are the steady-state density in the molecular region and the steady-state longitudinal current. It is therefore possible to obtain both these quantities by a steady-state self-consistent calculation without any need of knowing how the system has attained the steady state. i-DFT reduces to the LB+DFT approach if the xc bias correction is discarded and if the Hxc potential is the one of standard DFT depending only on the density.

Of course, for i-DFT to prove useful for the calculation of the transport properties of strongly correlated junctions one needs accurate approximations to the i-DFT functionals. Unfortunately, we are still not in the position of offering (H)xc potentials ready to use in available first-principles codes. Yet, we have been able to identify crucial and general properties that any good approximation to these functionals should fulfill. We have shown that the step of the equilibrium Hxc potential as the total number of particles crosses an integer bifurcates at finite current. Furthermore, a bifurcating step structure at integer number of particles occurs also for the xc bias as the current starts flowing. This latter property is of utmost importance since the derivative $\partial V_{\rm xc}/\partial I$ enters explicitly into the i-DFT formula for the zero-bias conductance and it plays a crucial role in suppressing the Kondo plateau in the CB regime. We have also identified an interesting duality between the values of the current and number of particles inside a CB diamond on one side and the intersections between discontinuity lines in the (H)xc potentials on the other side. These qualitative features should be incorporated in any approximation in order to reproduce the CB pattern of the differential conductance.

As a proof of concept we have examined in detail the SIAM and provided an accurate parametrization of the i-DFT potentials by best fitting the numerical results from fRG and NRG. The performance of i-DFT has turned out to be a full success. The differential conductance at any finite bias calculated by solving the i-DFT equations is accurate for any value of the interaction strenght and for temperatures ranging from zero to well above the Kondo temperature. Last but not least, owing to the simplicity of the i-DFT equations the calculation of an $I-V$ characteristics takes less than a CPU second.

An interesting possible extension of i-DFT which we are currently exploring consists in considering the Hxc potential and xc bias as adiabatic functionals of current and density for time-dependent calculations. This study generalizes previous investigations in strongly correlated models where the equilibrium xc potential of DFT LimaOliveiraCapelle:02 ; LimaSilvaOliveiraCapelle:03 ; SchoenhammerGunnarssonNoack:95 was turned, through the adiabatic approximation, into the xc potential of TDDFT Verdozzi:08 ; kskvg.2010 ; UimonenKhosraviStanStefanucciKurthLeeuwenGross:11 ; KhosraviUimonenStanStefanucciKurthLeeuwenGross:12 . In Ref. UimonenKhosraviStanStefanucciKurthLeeuwenGross:11, it was pointed out that the inaccurate value of the TDDFT current in the SIAM was due to the neglect of the xc bias. However, as we have shown, this is not the whole story. It is crucial that the xc bias is also a functional of the current. In fact, the time-local dependence on the current in i-DFT translates into a time- and space-nonlocal dependence on the density in TDDFT, see Ref. nptvc.2007, . The entangled space and time nonlocality of the functionals are strongly related by conservation laws vk.1996 . Thus, i-DFT holds promise for an improved description of time-dependent phenomena like, e.g., transient processes or AC responses.

Acknowledgements.

S.K. acknowledges funding by a grant of the ”Ministerio de Economia y Competividad (MINECO)” (FIS2016-79464-P) and by the “Grupos Consolidados UPV/EHU del Gobierno Vasco” (IT578-13). G.S. acknowledges funding by MIUR FIRB Grant No. RBFR12SW0J and EC funding through the RISE Co-ExAN (GA644076).

Appendix A Exact equilibrium density of the Constant Interaction Model

In this Appendix we describe an algorithm which allows for the calculation of the equilibrium density of the (uncontacted) Constant Interaction Model (CIM) described by the Hamiltonian $\hat{H}^{\rm CIM}$ of Eq. (40) at arbitrary temperature. Our aim is to calculate the densities (occupations)

[TABLE]

where

[TABLE]

with the partition function $Z_{K}(\tilde{\varepsilon}_{1},\ldots,\tilde{\varepsilon}_{K})$ of the CIM with $K$ single-particle levels with energies $\tilde{\varepsilon_{k}}=\varepsilon_{k}-\mu$. Suppose that we know the partition function $Z_{K-1}(\tilde{\varepsilon}_{1},\ldots,\tilde{\varepsilon}_{K-1})$ of the system with $K-1$ levels. Then the partition function for $K$ levels can be calculated by the recursion relation

[TABLE]

Defining the quantity

[TABLE]

from Eq. (147) one can then easily derive the recursive relation

[TABLE]

where we have also defined

[TABLE]

With these definitions the occupation of level $k$ becomes

[TABLE]