Non-equilibrium spectral functions from multi-terminal steady-state density functional theory

TL;DR

This paper extends steady-state density functional theory to multi-terminal nanoscale systems, enabling the extraction of non-equilibrium spectral functions and capturing phenomena like Kondo resonance splitting under bias.

Contribution

It develops a multi-terminal i-DFT formalism that can describe complex transport setups and non-equilibrium spectral features in nanoscale junctions.

Findings

Successfully extracts non-equilibrium spectral functions in a three-terminal setup.

Capable of describing Kondo resonance splitting under bias.

Validates the formalism against numerically exact many-body approaches.

Abstract

Multi-terminal transport setups allow to realize more complex measurements and functionalities (e.g., transistors) of nanoscale systems than the simple two-terminal arrangement. Here the steady-state density functional formalism (i-DFT) for the description of transport through nanoscale junctions with an arbitrary number of leads is developed. In a three-terminal setup and in the ideal STM limit where one of the electrodes (the `STM tip') is effectively decoupled from the junction, the formalism allows to extract its non-equilibrium spectral function (at arbitrary temperature) while a bias is applied between the other two electrodes. Multi-terminal i-DFT is shown to be capable of describing the splitting of the Kondo resonance in an Anderson impurity in the presence of an applied bias voltage, as predicted by numerically exact many-body approaches.

Equations

$$I_{\mathrm{T}}=2\int \mathrm{Tr}\left[f_{\mathrm{T}}(\omega -V_{\mathrm{T}}){\Gamma }_{\mathrm{T}}A(\omega )+i{\Gamma }_{\mathrm{T}}{G}^{<}(\omega )\right]$$

$$\underset{{\Gamma }_{\mathrm{T}}\to 0}{\lim }\frac{\partial G^{\lessgtr }(\omega )}{\partial V_{\mathrm{T}}}=0.$$

$$\Gamma_{\rm T}=\gamma_{{\rm T}}\bigg{[}\eta_{p}|p\rangle\langle p|+\eta_{q}|q\rangle\langle q|+\eta_{pq}\big{(}|p\rangle\langle q|+|q\rangle\langle p|\big{)}\bigg{]}$$

$$\lim_{\gamma_{{\rm T}}\to 0}\frac{1}{\gamma_{{\rm T}}}\frac{\partial I_{{\rm T}}}{\partial V_{{\rm T}}}=\frac{\mbox{$\mathcal{A}$}(V_{{\rm T}})}{\pi}$$

$$\mbox{$\mathcal{A}$}(\omega)=\eta_{p}A_{pp}(\omega)+\eta_{q}A_{qq}(\omega)+\eta_{pq}\big{[}A_{pq}(\omega)+A_{qp}(\omega)\big{]}$$

$$n(\mathbf{r})$$

$$I_{\alpha }$$

$$\underset{{\Gamma }_{\mathrm{T}}\to 0}{\lim }\frac{\partial v_{\mathrm{Hxc}}}{\partial V_{\mathrm{T}}}=\underset{{\Gamma }_{\mathrm{T}}\to 0}{\lim }\frac{\partial V_{\mathrm{xc}}}{\partial V_{\mathrm{T}}}=0.$$

$$\lim_{\gamma_{{\rm T}}\to 0}\frac{1}{\gamma_{{\rm T}}}\frac{\partial I_{{\rm T}}}{\partial V_{{\rm T}}}=\frac{\mbox{$\mathcal{A}$}_{s}(V_{{\rm T}}+V_{{\rm T},\rm xc})}{\pi}\left(1+\frac{\partial V_{{\rm T},\rm xc}}{\partial I_{{\rm T}}}\frac{\partial I_{{\rm T}}}{\partial V_{{\rm T}}}\right)$$

$$\mbox{$\mathcal{A}$}(\omega)=\lim_{\gamma_{{\rm T}}\to 0}\frac{\mbox{$\mathcal{A}$}_{s}(\omega+V_{{\rm T},\rm xc}(\omega))}{1-\frac{\gamma_{{\rm T}}}{\pi}\frac{\partial V_{{\rm T},\rm xc}(\omega)}{\partial I_{{\rm T}}}\mbox{$\mathcal{A}$}_{s}(\omega+V_{{\rm T},\rm xc}(\omega))},$$

$$N=2\int \sum _{\alpha =\mathrm{L},\mathrm{R},\mathrm{T}}{\tilde{f}}_{\alpha }(\omega )\frac{{\gamma }_{\alpha }}{\gamma }{A}_s(\omega )$$

$$I_{\mathrm{T}}=2{\gamma }_{\mathrm{T}}\int \left[\frac{{\gamma }_{\mathrm{L}}+{\gamma }_{\mathrm{R}}}{\gamma }{\tilde{f}}_{\mathrm{T}}(\omega )-\sum _{\alpha =\mathrm{L},\mathrm{R}}\frac{{\gamma }_{\alpha }}{\gamma }{\tilde{f}}_{\alpha }(\omega )\right]A_s(\omega )$$

$$I=2\int [{\gamma }_{\mathrm{L}}\frac{2{\gamma }_{\mathrm{R}}+{\gamma }_{\mathrm{T}}}{\gamma }{\tilde{f}}_{\mathrm{L}}(\omega )-{\gamma }_{\mathrm{R}}\frac{2{\gamma }_{\mathrm{L}}+{\gamma }_{\mathrm{T}}}{\gamma }{\tilde{f}}_{\mathrm{R}}(\omega )$$

$$+\frac{{\gamma }_{\mathrm{T}}({\gamma }_{\mathrm{L}}-{\gamma }_{\mathrm{R}})}{\gamma }{\tilde{f}}_{\mathrm{T}}(\omega )]A_s(\omega )$$

$$N+\frac{I_{\alpha }}{{\gamma }_{\alpha }}=2\int f_{\alpha }(\omega -V_{\alpha })A(\omega ).$$

$$A(\omega )=\frac{N}{2}{l}_{\gamma }(\omega -v-U)+\left(1-\frac{N}{2}\right)l_{\gamma }(\omega -v).$$

$$v_{\mathrm{Hxc}}^{\mathrm{CB}}-V_{\alpha ,\mathrm{xc}}^{\mathrm{CB}}\approx \frac{U}{2}+\frac{U}{\pi }\mathrm{atan}\left[\frac{N+I_{\alpha }/(2{\gamma }_{\alpha })-1}{\nu W({\Theta }_{\alpha })}\right]$$

$$\displaystyle\mbox{$\mathcal{G}$}_{\alpha\beta}=\frac{\partial I_{\alpha}}{\partial V_{s,\mu}}\frac{\partial V_{s,\mu}}{\partial V_{\beta}}=\mbox{$\mathcal{G}$}_{s,\alpha\mu}\left(\delta_{\mu\beta}+\frac{\partial V_{\mu,\rm xc}}{\partial I_{\nu}}\mbox{$\mathcal{G}$}_{\nu\beta}\right)$$

$$V_{{\rm T},\rm xc}(N,I_{{\rm T}},I)=\big{[}1-b(N)a_{{\rm T}}(I_{{\rm T}})a(I)\big{]}V^{\rm CB}_{{\rm T},\rm xc}(N,I_{{\rm T}},I)$$

$$a_{\mathrm{T}}(I_{\mathrm{T}})=1-\frac{2}{\pi }\mathrm{atan}\left[\lambda {\left(\frac{I_{\mathrm{T}}}{W(0){\gamma }_{\mathrm{T},\mathrm{eff}}}\right)}^2\right]$$

$$a(I)=1-\frac{2}{\pi }\mathrm{atan}\left[\lambda {\left(\frac{I}{W(\Theta ){\gamma }_{\mathrm{eff}}}\right)}^2\right]$$

$$b(N=1)=1+\frac{1}{\left.\frac{\partial V^{\rm CB}_{{\rm T},\rm xc}}{\partial I_{{\rm T}}}\right|_{\begin{subarray}{c}N=1\\ I=I_{{\rm T}}=0\end{subarray}}}\left(\frac{1}{\mbox{$\mathcal{G}$}_{{\rm T}}}-\frac{1}{\mbox{$\mathcal{G}$}_{s,{\rm T}}}\right).$$

$$A(\omega )=\frac{\mathrm{d}}{\mathrm{d}\omega }{\int }^{\omega +V_{\mathrm{T},\mathrm{xc}}(\omega )}\mathrm{d}{\omega }^{\prime }{A}_s({\omega }^{\prime }).$$

Full text

I Introduction

Nanoscale or molecular junctions can now be made in the lab by connecting, e.g., individual atoms, molecules or clusters to electrodes, which may become the building blocks for prospective applications in Molecular Electronics Cuniberti et al. (2005); Cuevas and Scheer (2010) and/or Quantum TechnologiesJ.P.Dowling and G.J.Milburn (2003). On a more fundamental level, nanoscale junctions under an applied bias voltage are used to experimentally Cuevas and Scheer (2010); Ward et al. (2011); Tewari and van Ruitenbeek (2018) study the many-body problem of interacting electrons driven out of equilibrium, and to probe nanoscale quantum systems by differential conductance $dI/dV$ spectroscopy Wiesendanger (1994); Hamers and Padowitz (2001); Madhavan et al. (1998); Hirjibehedin et al. (2007).

The description of transport through real nanoscale junctions presents a major theoretical challenge since quantum effects, the atomistic details of the junction, electronic interactions and the out-of-equilibrium situation have to be properly taken into account. In its simplest form a nanoscale junction can be described by a single impurity Anderson model (SIAM) coupled to two leads at different chemical potentials that define the bias voltage across the junction. An intriguing prediction for this model is the splitting of the Kondo resonance under an applied bias voltage Meir et al. (1993); Wingreen and Meir (1994); Sun and Guo (2001); Krawiec and Wysokiński (2002); Shah and Rosch (2006); Fritsch and Kehrein (2010); Cohen et al. (2014). Experimentally, this effect cannot be seen directly in the $dI/dV$ of a two terminal device, as it only shows up in the non-equilibrium spectral function of the junction. However, it can be measured in a three-terminal setup where one of the electrodes is very weakly coupled and serves as a probe Lebanon and Schiller (2001); Sun and Guo (2001); De Franceschi et al. (2002); Leturcq et al. (2005).

In general, many-body techniques for solving the out-of-equilibrium problem are computationally too demanding to be applied to more than relatively simple model systems such as the SIAM and slightly more complex models. Owing to its conceptual simplicity and computational efficiency, the now standard approach for realistic modeling of electronic transport in nanoscale junctions combines density functional theory (DFT) calculations with the Landauer-Büttiker approach (LB) to transport Taylor et al. (2001); Palacios et al. (2002). While the LB-DFT approach properly takes into account atomistic details of the junctions as well as quantum effects, it is formally incomplete in the sense that there is no guarantee that it gives the correct current through the interacting system even if the exact exchange-correlation functional is used Kurth and Stefanucci (2013); Stefanucci and Kurth (2013). It is therefore not surprising, that the LB-DFT formalism does not capture all aspects of correlated electronic transport, namely Coulomb blockade and Kondo physics, although under special circumstances some of these aspects may be correctly described in a surprisingly simple manner Stefanucci and Kurth (2011); Bergfield et al. (2012); Tröster et al. (2012). Combination of the LB-DFT approach with many-body methods incorporates electronic correlations (originating from a relatively small subspace) into the description of electronic transport through realistic systems Jacob et al. (2009); Jacob and Kotliar (2010); Jacob (2015); Droghetti and Rungger (2017), but suffers from the infamous double-counting problem.

Recently, a novel approach, called steady-state DFT (or i-DFT), has been devised to describe the steady-state transport through nanoscale junctions driven out of equilibrium in a DFT framework Stefanucci and Kurth (2015); Kurth and Stefanucci (2016, 2017). In i-DFT the steady current through the nanoscale junction is an additional fundamental “density” variable and the bias voltage across the junction is the corresponding potential. Provided that good approximations for the functionals are found, this approach is able to describe the full phenomenology of electronic correlations in transport through nanoscale junctions. Moreover, the i-DFT formalism can be applied to extract the equilibrium many-body spectral function from a DFT calulation Jacob and Kurth (2018).

In this Letter we generalize i-DFT to the multi-terminal situation, and then consider the specific situation of a junction connected to three electrodes. We show how in the “ideal STM setup” where one of the electrodes is only weakly coupled to the system, one can extract the non-equilibrium many-body spectral function of the junction at arbitrary temperature and bias between the other two electrodes. We apply the approach to the SIAM for which we construct an approximate xc functional which partially captures the splitting of the Kondo peak under finite bias. We also identify the crucial feature of the xc functional needed to fully describe the splitting of the Kondo peak.

II Non-equilibrium spectral functions

Here we breifly recall how to calculate out-of-equilibrium spectral functions from transport measurements Meir et al. (1993); Wingreen and Meir (1994); Sun and Guo (2001); Krawiec and Wysokiński (2002); Shah and Rosch (2006); Fritsch and Kehrein (2010); Cohen et al. (2014). We consider a three-terminal molecular junction as illustrated in Fig. 1. Two electrodes, the left ($\alpha={\rm L}$) and right ($\alpha={\rm R}$) ones, have voltages $V_{{\rm L}}=-V_{{\rm R}}=V/2$ (gauge fixing) and the same finite temperature $\Theta_{{\rm L}}=\Theta_{{\rm R}}=\Theta$. The third electrode plays the role of a tip ($\alpha={\rm T}$) and is kept at zero temperature and voltage $V_{{\rm T}}$. The contact between the tip and the nanoscopic region is described by the energy-independent hybridization $\Gamma_{{\rm T}}$ whose indices run over a suitable one-electron orbital basis for the considered molecule. The $\Gamma_{{\rm T}}$ matrix, aside from being constrained to be symmetric and positive semi-definite, will be varied at will.

According to Meir and Wingreen Meir and Wingreen (1992) the current $I_{{\rm T}}$ flowing out of the tip is given by (henceforth $\int\equiv\int\frac{d\omega}{2\pi}$)

[TABLE]

where $A(\omega)=i\left[G^{>}(\omega)-G^{<}(\omega)\right]$ is the nonequilibrium, finite-temperature many-body spectral function in terms of the lesser/greater Green’s functions whereas $f_{{\rm T}}(\omega)=\theta(-\omega)$ is the zero-temperature Fermi function of the tip. In Eq. (1) the trace is over the indices of the molecular one-electron basis. Similar to what we showed in previous work Jacob and Kurth (2018) in the ideal Scanning Tunneling Microscopy (STM) limit, $\Gamma_{{\rm T}}\to 0$, the Green’s functions $G^{\lessgtr}$ are not affected by a change of the tip voltage and hence

[TABLE]

We then consider a hybridization of the form

[TABLE]

This operator is symmetric and positive semi-definite for all $\gamma_{\rm T},\eta_{p},\eta_{q}>0$ and $|\eta_{pq}|\leq\sqrt{\eta_{p}\eta_{q}}$. Taking into account Eq. (2) it is straightforward to show that

[TABLE]

where

[TABLE]

is a linear combination of the matrix elements of the spectral function, i.e. $A_{pq}(\omega)=\langle p|A(\omega)|q\rangle$. Choosing, e.g., $\eta_{p}=1$ and $\eta_{q}=0$ we can obtain all diagonal elements A_{pp}=\mbox{\mathcal{A}} by varying $p$. Subsequently we can extract the off-diagonal elements A_{pq}+A_{qp}=\mbox{\mathcal{A}}-A_{pp}-A_{qq} by setting $\eta_{pq}=\eta_{p}=\eta_{q}=1$.

III Multi-terminal i-DFT

In Ref. Jacob and Kurth (2018) we showed how to calculate equilibrium and zero-temperature spectral functions from the i-DFT approach Stefanucci and Kurth (2015). For nonequilibrium and finite-temperature spectral functions we have to generalize i-DFT to multi-terminal setups, with electrodes at different voltages and temperatures.

We consider a nanoscopic region $\mathcal{M}$ containing a quantum dot or molecule and a number of electrodes \alpha=1,..,\mbox{\mathcal{N}}, as depicted schematically in Fig. 1 for \mbox{\mathcal{N}}=3. The system is assumed to be in a steady state characterized by temperatures $\Theta_{\alpha}$ and external voltages $V_{\alpha}$ in electrode $\alpha$ and by a gate voltage $v({\mathbf{r}})$ in $\mathcal{M}$. As long as region $\mathcal{M}$ is finite there are no constraints on the shape of its boundaries. Due to gauge invariance the same steady-state is attained by shifting all voltages by a constant energy $W$, i.e., $V_{\alpha}\to V_{\alpha}+W$ and $v({\mathbf{r}})\to v({\mathbf{r}})+W$. Let $I_{\alpha}$ be the longitudinal current flowing out of electrode $\alpha$ and $n({\mathbf{r}})$ be the density in the nanoscopic region. Due to charge conservation (consequence of the aforementioned gauge invariance) the currents fulfill $\sum_{\alpha}I_{\alpha}=0$. With a similar proof as the one published in Ref. Stefanucci and Kurth (2015) we can state the multi-terminal generalization of the i-DFT theorem:

Theorem: There exists a one-to-one map between the set of “densities” (n,I_{1},\ldots,I_{\tiny\mbox{\mathcal{N}}}) with $\sum_{\alpha}I_{\alpha}=0$ and the set of “potentials” (v,V_{1},\ldots,V_{\tiny\mbox{\mathcal{N}}}) up to a constant shift $W$. The bijectivity of the map is guaranteed in a finite (and gate dependent) region around zero voltages $V_{\alpha}$ for any set of finite temperatures $\Theta_{\alpha}$.

According to the multi-terminal i-DFT theorem there exists a unique set of Kohn-Sham (KS) potentials $(v_{s}$, V_{1,s},\ldots,V_{\tiny\mbox{\mathcal{N}},s}) which in the noninteracting system reproduce the density $n({\mathbf{r}})$ and currents \mbox{\mathcal{I}}=(I_{1},\ldots,I_{\tiny\mbox{\mathcal{N}}}) of the interacting system (here we are assuming that the density and the currents are non-interacting representable). Following the KS procedure we define the exchange-correlation (xc) voltages V_{\alpha,\rm xc}[n,\mbox{\mathcal{I}}]=V_{\alpha,s}[n,\mbox{\mathcal{I}}]-V_{\alpha}[n,\mbox{\mathcal{I}}] and the Hartree-xc (Hxc) gate voltage v_{\rm Hxc}[n,\mbox{\mathcal{I}}]=v_{s}[n,\mbox{\mathcal{I}}]-v[n,\mbox{\mathcal{I}}] (which are functionals of the density in $\mathcal{M}$ and the currents) and then calculate the interacting density and currents by solving self-consistently the equations

[TABLE]

where $f_{\alpha}(\omega)=1/(e^{\omega/\Theta_{\alpha}}+1)$ is the Fermi function of lead $\alpha$ at temperature $\Theta_{\alpha}$. In the KS equations $A_{\alpha,s}(\omega,{\mathbf{r}})=\langle{\mathbf{r}}|G^{\rm R}_{s}(\omega)\Gamma_{\alpha}(\omega)G^{\rm A}_{s}(\omega)|{\mathbf{r}}\rangle$ is the partial KS spectral function written in terms of the retarded/advanced KS Green’s functions $G^{\rm R/A}_{s}$ and hybridization $\Gamma_{\alpha}(\omega)$ due to lead $\alpha$, whereas \mbox{\mathcal{T}}_{\alpha\alpha^{\prime},s}(\omega)={\rm Tr}\left[G^{\rm R}_{s}(\omega)\Gamma_{\alpha}(\omega)G^{\rm A}_{s}(\omega)\Gamma_{\alpha^{\prime}}(\omega)\right] are the KS transmission probabilities.

IV Spectral function from i-DFT

We specialize the multi-terminal i-DFT theorem to the three-terminal case previously discussed. Let us fix the gauge according to $V_{{\rm L},s}=-V_{{\rm R},s}=V_{s}/2$ and let us consider the combination $I=(I_{L}-I_{R})/2$ and $I_{{\rm T}}$ as the two independent currents. Then the triple $v_{\rm Hxc}=v_{\rm Hxc}[n,I_{{\rm T}},I]$, $V_{{\rm T},\rm xc}=V_{{\rm T},\rm xc}[n,I_{{\rm T}},I]$ and $V_{\rm xc}=V_{\rm xc}[n,I_{{\rm T}},I]$ are functionals of the triple $n$, $I_{{\rm T}}$ and $I$ (here $V_{\rm xc}[n,I_{{\rm T}},I]=V_{s}[n,I_{{\rm T}},I]-V[n,I_{{\rm T}},I]$). Considering $n,I_{{\rm T}}$ and $I$ as interacting functionals of the physical voltages $v,V_{{\rm T}}$ and $V$, Eq. (2) implies that $\partial{n({\mathbf{r}})}/\partial{V_{\rm T}}\rightarrow 0$ and $\partial{I}/\partial{V_{\rm T}}\rightarrow 0$ for $\Gamma_{\rm T}\rightarrow 0$, and by the chain rule it thus follows that

[TABLE]

In the same spirit as in our previous work Jacob and Kurth (2018), we now take advantage of these relations in order to express the spectral function $A$ in terms of the KS spectral function $A_{s}$. In the noninteracting KS system the tip current is given by Eq. (1), replacing $A(\omega)$ with the KS spectral function $A_{s}=\sum_{\alpha}A_{\alpha,s}$ and $G^{<}(\omega)$ by the KS lesser GF $G^{<}_{s}$. Taking into account Eq. (8) and the fact that $\Gamma_{{\rm T}}$ is energy-independent we find

[TABLE]

where \mbox{\mathcal{A}}_{s} is defined as in Eq. (5) with $A\to A_{s}$. Combining this result with Eq. (4) we arrive at the first main result of this work

[TABLE]

which generalizes the corresponding result of Ref. Jacob and Kurth (2018) to nonequilibrium spectral functions. Here we have made explicit the dependence of $V_{{\rm T},\rm xc}$ on $\omega=V_{{\rm T}}$ through its dependence on $I_{{\rm T}}$. Choosing, e.g., $\eta_{p}=1$ and $\eta_{q}=0$, Eq. (10) provides a relation between $A_{pp}$ and $A_{s,pp}$. The off-diagonal combination $A_{pq}+A_{qp}$ does instead follow by setting $\eta_{pq}=\eta_{p}=\eta_{q}=1$. We also observe that both \mbox{\mathcal{A}}_{s} and $\mathcal{A}$ are normalized to the same value, i.e. \int\mbox{\mathcal{A}}(\omega)=\int\mbox{\mathcal{A}}_{s}(\omega) as it should be 111 This follows by integrating over $\omega$ both sides of Eq. (10), changing variable $\omega^{\prime}=\omega+V_{{\rm T},\rm xc}(\omega)$ in the r.h.s. and taking into account the Jacobian \frac{{\rm d}\omega^{\prime}}{{\rm d}\omega}=1/(1-\frac{\gamma_{{\rm T}}}{\pi}\frac{\partial V_{\rm T,xc}(\omega)}{\partial I_{T}}\mbox{\mathcal{A}}_{s}(\omega+V_{\rm T,xc}(\omega))). .

V i-DFT potentials for the Anderson model

We apply the i-DFT framework to the single-impurity Anderson model (SIAM) with charging energy $U$. Since the SIAM nanoscopic region has only one electronic degree of freedom the density $n=N$ coincides with the impurity occupation $N$, and all hybridization matrices are scalar. We then write $\Gamma_{\rm T}=\gamma_{{\rm T}}$ for the tip and consider energy-independent left/right hybridizations $\gamma_{{\rm L}/{\rm R}}$. The i-DFT self-consistent equations for $N$, $I_{{\rm T}}$ and $I$ read

[TABLE]

[TABLE]

[TABLE]

where we have defined $\tilde{f}_{\alpha}(\omega)\equiv f_{\alpha}(\omega-V_{\alpha}-V_{\alpha,\rm xc})$ as the shifted Fermi function and $\gamma\equiv\gamma_{{\rm L}}+\gamma_{{\rm R}}+\gamma_{{\rm T}}$. The KS spectral function is simply $A_{s}(\omega)=\ell_{\gamma}(\omega-v-v_{\rm Hxc})$ with the Lorentzian $\ell_{\gamma}(\omega)=\gamma/(\omega^{2}+\gamma^{2}/4)$.

In order to derive an approximation for the i-DFT potentials we observe that in the interacting system the current flowing out of lead $\alpha$ reads $I_{\alpha}=2\int\left[f_{\alpha}(\omega-V_{\alpha})\gamma_{\alpha}A(\omega)+i\gamma_{\alpha}G^{<}(\omega)\right]$. Taking into account that the impurity occupation is $N=-2i\int\frac{d\omega}{2\pi}G^{<}(\omega)$ we get

[TABLE]

In the CB regime, i.e., for temperatures $\Theta_{{\rm L}}=\Theta_{{\rm R}}=\Theta$ larger than the Kondo temperature (at ph symmetry Jakobs et al. (2010)) $\Theta_{K}=\frac{4}{\pi}\sqrt{U\gamma}\exp\left(-\frac{\pi}{4}\left(\frac{U}{\gamma}-\frac{\gamma}{U}\right)\right)$ but smaller than $\gamma_{{\rm L}/{\rm R}}$, the interacting spectral function is well approximated by Stefanucci and Kurth (2015); Dittmann et al. (2018)

[TABLE]

Inserting Eq. (15) in the r.h.s. of Eq. (14) we get the same expression obtained in Ref. Stefanucci and Kurth (2015) for the two-terminal set-up. Therefore, the CB reverse-engineered xc potentials can be parametrized in the same manner

[TABLE]

with $\nu=1$, $W(\Theta)=0.16\times(\gamma/U)(1+9(\Theta/\gamma)^{2})$ and $I_{{\rm L}}=I-I_{{\rm T}}/2$, $I_{{\rm R}}=-I-I_{{\rm T}}/2$ (as follows from charge conservation). From Eqs. (16) we can easily extract an explicit form of the (H)xc potentials $v^{\rm CB}_{\rm Hxc}$, $V^{\rm CB}_{{\rm T},\rm xc}$ and $V^{\rm CB}_{\rm xc}=2V^{\rm CB}_{{\rm L},\rm xc}=-2V^{\rm CB}_{{\rm R},\rm xc}$ in terms of $N,I_{{\rm T}}$ and $I$.

The (H)xc potentials in Eq. (16) are certainly inadequate for temperatures $\Theta\lesssim\Theta_{K}$. In particular for $\Theta=0$ the Friedel sum rule implies that the zero-bias interacting and KS conductances \mbox{\mathcal{G}}_{\alpha\beta}=\partial I_{\alpha}/\partial V_{\beta} and \mbox{\mathcal{G}}_{s,\alpha\beta}=\partial I_{\alpha}/\partial V_{s,\beta} are identical 222Using the Friedel sum-rule one can show that \mbox{\mathcal{G}}_{\alpha\beta}=C_{\alpha\beta}(\gamma_{\rm T},\gamma_{\rm L},\gamma_{\rm R})A(\mu) where the prefactor $C_{\alpha\beta}$ depends only on the hybridizations and $A(\mu)$ is the interacting spectral function at chemical potential $\mu$ (which is set to zero in our case). Since $A(\mu)=\frac{4}{\gamma}\sin^{2}(\pi N/2)$ and since in i-DFT the KS occupation $N$ is the same as the interacting $N$ we conclude that the interacting and KS conductances are the same.. Since (repeated indices are summed over)

[TABLE]

the zero-temperature xc voltages must fulfill $\partial V_{\mu,\rm xc}/\partial I_{\nu}=0$ at zero currents. We incorporate this property in $v_{\rm Hxc}$ and $V_{\rm xc}$ using the same parametrization proposed in Ref. Stefanucci and Kurth (2015) for the two-terminal case, i.e., for $\gamma_{{\rm T}}=0$, which has been shown to be accurate in a wide range of temperatures and charging energy. For $V_{{\rm T},\rm xc}$ we propose

[TABLE]

where in $V^{\rm CB}_{{\rm T},\rm xc}$ we now take $\nu=2$ Kurth and Stefanucci (2016) and the functions $a_{{\rm T}}$ and $a$ are similar to the one used in Ref. Jacob and Kurth (2018) and read

[TABLE]

[TABLE]

with $\gamma_{{\rm T},\rm eff}=\frac{4\gamma_{{\rm T}}(\gamma_{{\rm L}}+\gamma_{{\rm R}})}{\gamma}$, $\gamma_{\rm eff}=\frac{4\gamma_{{\rm L}}\gamma_{{\rm R}}}{\gamma_{{\rm L}}+\gamma_{{\rm R}}}$ and $\lambda=0.16$ the same fit parameter used in Ref. Jacob and Kurth (2018). For $b(N)$ we implement the same function as in Ref. Stefanucci and Kurth (2015) but we replace the two-terminal conductance \mbox{\mathcal{G}}_{\rm univ}=dI/dV at the ph symmetric gate $v=-U/2$, voltage $V=0$ and symmetric coupling $\gamma_{{\rm L}}=\gamma_{{\rm R}}$ (this is a universal function depending only on the ratio $\Theta/\Theta_{K}$) with the three-terminal conductance \mbox{\mathcal{G}}_{{\rm T}}=dI_{{\rm T}}/dV_{{\rm T}} at the ph symmetric gate and voltages $V=V_{{\rm T}}=0$:

[TABLE]

One can show that \mbox{\mathcal{G}}_{{\rm T}}=\frac{4\gamma_{{\rm T}}(\gamma_{{\rm L}}+\gamma_{{\rm R}})}{\gamma^{2}}\mbox{\mathcal{G}}_{\rm univ}. In Eq. (21) \mbox{\mathcal{G}}_{s,{\rm T}}=dI_{{\rm T}}/dV_{s,{\rm T}} is the KS conductance at the same external potentials, i.e., ph gate and zero voltages.

VI Results

As a first test we use our three-terminal i-DFT setup to compute the spectral function of the SIAM in thermal equilibrium for which we can compare with results from numerical renormalization group (NRG) techniques Motahari (2017); Requist , see Fig. 2. The i-DFT spectra agree reasonably well with the NRG ones although the height of the Kondo peak is slightly overestimated and for $\Theta/\Theta_{K}\gtrsim 2.5$ the Coulomb blockade side peaks are a bit too narrow. In general, the finite temperature i-DFT spectra are of comparable quality as the zero-temperature ones Jacob and Kurth (2018).

We now consider the zero-temperature, non-equilibrium SIAM and benchmark the i-DFT spectra against recent results from the Quantum Monte Carlo (QMC) approach Bertrand et al. (2019), see Fig. 3. i-DFT reproduces all main qualitative features of the QMC spectra. In particular, our simple functional of Eq. (18) for the xc tip bias is able to capture the finite-bias splitting of the Kondo peak in this moderately correlated case $U/\gamma=2.5$. Nevertheless, in i-DFT the splitting appears at somewhat higher biases and the distance between the peaks increases with bias faster than in QMC. We have done calculations for the same set of biases but at a finite temperature $\Theta/\Theta_{K}=0.6$ and observed no dramatic changes except for the suppression of the Kondo peak already at zero voltage.

In Fig. 4 (left panel) we compare i-DFT with QMC non-equilibrium spectral functions Bertrand et al. (2019) for a stronger interaction strength $U/\gamma=4$. Clearly our approximation to $V_{\rm T,xc}$ is missing a crucial feature since the Kondo splitting is totally absent in i-DFT. Below we highlight an exact property that $V_{\rm T,xc}$ must fulfill in order to capture the finite-bias splitting. The interacting spectral function in Eq. (10) can also be written as

[TABLE]

Therefore, given a many-body (e.g., QMC) spectral function $A(\omega)$, by integration of Eq. (22) one can reverse-engineer the xc tip bias $V_{\rm T,xc}$ which corresponds to the given $A$. In the upper right panel of Fig. 3 we extracted $V_{\rm T,xc}$ as function of $I_{\rm T}$ (for fixed values of $N$ and $I$) corresponding to the QMC spectral functions of the left panel of the same figure and compare to our i-DFT functional of Eq. (18). Although some differences are visible our approximate xc tip bias seems to agree rather well with the reverse engineered one. The missing feature becomes evident if we compare the derivatives of $V_{\rm T,xc}$ w.r.t. $I_{\rm T}$, see lower right panel of Fig. 3. While the derivative of the reverse engineered $V_{\rm T,xc}$ exhibits a double peak in the vicinity of $I_{\rm T}/\gamma_{T}\approx 0$, our approximation exhibits only a single maximum at $I_{\rm T}/\gamma_{T}=0$. Of course, the height as well as the positions of the maxima depend on the current $I$ between the left and right leads. We have verified that using the reverse engineered $V_{\rm T,xc}$ in Eq. (10) the i-DFT and QMC spectral functions become indistinguishable. The correct incorporation of the double peak feature into an improved approximation for $V_{\rm T,xc}$ is beyond the scope of this work. However, the established existence of this xc bias constitues a proof-of-concept: i-DFT provides a numerically cheap method to calculate non-equilibrium spectral functions at zero and finite temperature.

VII Conclusions

We have generalized the i-DFT formalism for steady state transport through nanoscale junctions to the situation of multiple electrodes. In particular, for a three-terminal setup in the limit of vanishing coupling to one of the electrodes (ideal STM limit) we have shown how to extract the non-equilibrium spectral function of the junction at both zero and finite temperature extending earlier work Jacob and Kurth (2018) which was restricted both to equilibrium and zero temperature. For the specific situation of an Anderson model coupled to three electrodes we have constructed an approximate xc functional by a relatively simple “natural” generalization of already existing i-DFT functionals. This approximation describes, at least for not too strong interactions, the splitting of the Kondo peak at finite bias and yields results in reasonable qualitative agreement with computationally more demanding many-body approaches such as NRG and non-equilibrium QMC. Although for stronger interactions our approximation does not capture the splitting of the Kondo peak, we were nevertheless able to identify the missing feature which needs to be incorporated in future improved functionals. In order to construct such functionals, reliable reference results from other many-body methods are certainly very welcome Cohen et al. (2014); Bertrand et al. (2019); Krivenko et al. (2019). However, once such approximations are available for a relatively simple system such as the Anderson model, generalizations to more complicated model systems (such as, e.g., multi-level systems) may actually be relatively straightforward Stefanucci and Kurth (2015); Kurth and Stefanucci (2017); Jacob and Kurth (2018); Kurth and Jacob (2018). Since multi-terminal i-DFT is comparable in computational effort to standard LB-DFT calculations, it is therefore suitable to study systems currently inaccessible for accurate out-of-equilibrium many-body methods.

S.K. acknowledges funding through a grant of the ”Ministerio de Economia y Competividad (MINECO)” (FIS2016-79464-P). GS acknowledges EC funding through the RISE Co-ExAN (Grant No. GA644076) and Tor Vergata University for financial support through the Mission Sustainability Project 2DUTOPI.