Absence of Coulomb Blockade in the Anderson Impurity Model at the Symmetric Point

TL;DR

This paper demonstrates that at the symmetric point, the Anderson impurity model exhibits no Coulomb blockade effects, with the steady-state current indistinguishable from a non-interacting model, challenging expectations of many-body interactions.

Contribution

It shows that at the symmetric point, the Anderson impurity model's steady-state current mimics a non-interacting system, revealing the absence of Coulomb blockade effects.

Findings

Steady-state current matches non-interacting models at the symmetric point.

Both methods accurately predict current in Coulomb blockade regime.

Many-body interactions do not influence the current at this symmetry point.

Abstract

In this work, we investigate the characteristics of the electric current in the so-called symmetric Anderson impurity model. We study the nonequilibrium model using two complementary approximate methods, the perturbative quantum master equation approach to the reduced density matrix, and a self-consistent equation of motion approach to the nonequilibrium Green's function. We find that at a particular symmetry point, an interacting Anderson impurity model recovers the same steady-state current as an equivalent non-interacting model, akin a two-band resonant level model. We show this in the Coulomb blockade regime for both high and low temperatures, where either the approximate master equation approach and the Green's function method provide accurate results for the current. We conclude that the steady-state current in the symmetric Anderson model at this regime does not encode characteristics of a many-body interacting system.

Figures

Equations

$$H_S=$$

$$H_B=$$

$$V=$$

$$\left(\begin{array}[]{c}\dot{p}_{0}\\ \dot{p}_{\uparrow}\\ \dot{p}_{\downarrow}\\ \dot{p}_{f}\end{array}\right)=M\left(\begin{array}[]{c}p_{0}\\ p_{\uparrow}\\ p_{\downarrow}\\ p_{f}\end{array}\right),$$

$$M=\sum_{\ell=\rm L,R}\left(\begin{array}[]{cccc}-\left(\gamma^{\ell}_{\uparrow 0}+\gamma^{\ell}_{\downarrow 0}\right)&\gamma^{\ell}_{0\uparrow}&\gamma^{\ell}_{0\downarrow}&0\\ \gamma^{\ell}_{\uparrow 0}&-\left(\gamma^{\ell}_{0\uparrow}+\gamma^{\ell}_{f\uparrow}\right)&0&\gamma^{\ell}_{\uparrow f}\\ \gamma^{\ell}_{\downarrow 0}&0&-\left(\gamma^{\ell}_{0\downarrow}+\gamma^{\ell}_{f\downarrow}\right)&\gamma_{\downarrow f}\\ 0&\gamma^{\ell}_{f\uparrow}&\gamma^{\ell}_{f\downarrow}&-\left(\gamma^{\ell}_{\uparrow f}+\gamma^{\ell}_{\downarrow f}\right)\end{array}\right),$$

$${\gamma }_{\uparrow 0}^{\ell }$$

$${\gamma }_{0\uparrow }^{\ell }$$

$${\gamma }_{f\uparrow }^{\ell }$$

$${\gamma }_{\uparrow f}^{\ell }$$

$$I_{\ell }(t)=\frac{e}{2}\sum _{ij}\left(N_i-N_j\right)J_{ij}^{\ell }\text{where}{J}_{ij}^{\ell }={\gamma }_{ij}^{\ell }{p}_j(t)-{\gamma }_{ji}^{\ell }{p}_i(t),$$

$$G_{\sigma }\left(\tau ,{\tau }^{\prime }\right)=-i⟨{\hat{T}}_c\left[d_{\sigma }\left(\tau \right)d_{\sigma }^{†}\left({\tau }^{\prime }\right)\right]⟩$$

$$G_{\sigma }^r\left(\omega \right)=$$

$$G_{\sigma }^{\left(2\right)r}\left(\omega \right)=$$

$$G_{\sigma }^{<}\left(\omega \right)=$$

$$+$$

$$G_{\sigma }^{\left(2\right)<}\left(\omega \right)=$$

$$+$$

$$G_{\sigma }^r\left(t,t^{\prime }\right)$$

$$G_{\sigma }^{<}\left(t,t^{\prime }\right)$$

$$G_{\sigma }^{\left(2\right)r}\left(t,t^{\prime }\right)$$

$$G_{\sigma }^{\left(2\right)<}\left(t,t^{\prime }\right)$$

$$⟨n_{\sigma }⟩={\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }{G}_{\sigma }^{<}\left(\omega \right),$$

$$g_{\sigma }^{\left(\delta \right)r}\left(\omega \right)$$

$$g_{\sigma }^{\left(\delta \right)<}\left(\omega \right)$$

$${\Sigma }_{\sigma }^r\left(\omega \right)$$

$${\Sigma }_{\sigma }^{<}\left(\omega \right)$$

$$I_{\ell }=2e\sum _{\sigma }\int \frac{d\omega }{2\pi }\left[i{\Sigma }_{\ell }^{<}\left(\omega \right)\mathrm{Im}\left\{G_{\sigma }^r\left(\omega \right)\right\}+{\Sigma }_{\ell }^r\left(\omega \right)G_{\sigma }^{<}\left(\omega \right)\right].$$

$$I_{\mathrm{L}}^{\left(U=0\right)}=\frac{2e{\Gamma }_{\mathrm{L}}{\Gamma }_{\mathrm{R}}}{{\Gamma }_{\mathrm{L}}+{\Gamma }_{\mathrm{R}}}\left(f^{\mathrm{L}}(\epsilon -\mu )-f^{\mathrm{R}}(\epsilon +\mu )\right),$$

$$I_{\mathrm{L}}^{\left(U=-2\epsilon \right)}=\frac{2e{\Gamma }_L{\Gamma }_R\left(f^{\mathrm{L}}(\epsilon -\mu )f^{\mathrm{R}}(\epsilon -\mu )-f^{\mathrm{L}}(\epsilon +\mu )f^{\mathrm{R}}(\epsilon +\mu )\right)}{{\Gamma }_{\mathrm{L}}\left(f^{\mathrm{L}}(\epsilon -\mu )+f^{\mathrm{L}}(\epsilon +\mu )\right)+{\Gamma }_{\mathrm{R}}\left(f^{\mathrm{R}}(\epsilon -\mu )+f^R(\epsilon +\mu )\right)}.$$

$$\displaystyle n_{\uparrow}\big{|}_{\begin{subarray}{c}U=0\end{subarray}}$$

$$\displaystyle(1-n_{\uparrow})\big{|}_{\begin{subarray}{c}U=0\end{subarray}}$$

$$I_{\mathrm{L}}^{\left(U=-2\epsilon \right)}-I_{\mathrm{L}}^{\left(U=0\right)}=e{\Gamma }^{2}U\frac{⟨n_{\uparrow }⟩+⟨n_{\downarrow }⟩}{2}{\int }_{-\infty }^{\infty }\frac{d\omega }{2\pi }\frac{2\omega \left(f^{\mathrm{L}}\left(\omega -{\mu }_{\mathrm{L}}\right)-f^{\mathrm{R}}\left(\omega -{\mu }_{\mathrm{R}}\right)\right)}{{\omega }^4+2{\omega }^2\left({\left(\frac{\Gamma }{2}\right)}^2-{\left(\frac{U}{2}\right)}^2\right)+{\left({\left(\frac{\Gamma }{2}\right)}^2+{\left(\frac{U}{2}\right)}^2\right)}^2}.$$

$$H_S=$$

$$H_B=$$

$$V=$$

$$\left(\begin{array}[]{c}\dot{p}_{0}\\ \dot{p}_{\uparrow}\\ \dot{p}_{\downarrow}\\ \dot{p}_{f}\end{array}\right)=M\left(\begin{array}[]{c}p_{0}\\ p_{\uparrow}\\ p_{\downarrow}\\ p_{f}\end{array}\right),$$

$$M=\sum_{\ell\in L,R}\left(\begin{array}[]{cccc}-\left(\gamma^{\ell}_{\uparrow 0}+\gamma^{\ell}_{\downarrow 0}\right)&\gamma^{\ell}_{0\uparrow}&\gamma^{\ell}_{0\downarrow}&0\\ \gamma^{\ell}_{\uparrow 0}&-\left(\gamma^{\ell}_{0\uparrow}+\gamma^{\ell}_{f\uparrow}\right)&0&\gamma^{\ell}_{\uparrow f}\\ \gamma^{\ell}_{\downarrow 0}&0&-\left(\gamma^{\ell}_{0\downarrow}+\gamma^{\ell}_{f\downarrow}\right)&\gamma_{\downarrow f}\\ 0&\gamma^{\ell}_{f\uparrow}&\gamma^{\ell}_{f\downarrow}&-\left(\gamma^{\ell}_{\uparrow f}+\gamma^{\ell}_{\downarrow f}\right)\end{array}\right).$$

Full text

\alsoaffiliation

The Raymond and Beverly Sackler Center for Computational Molecular and Materials Science, Tel Aviv University, Tel Aviv, Israel 69978 \altaffiliationBoth authors contributed equally to this work

\altaffiliationBoth authors contributed equally to this work

\alsoaffiliationKavli Energy NanoScience Institute, Berkeley, California 94720, United States \alsoaffiliationMaterials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, United States

\alsoaffiliationMaterials Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, United States \alsoaffiliationThe Raymond and Beverly Sackler Center for Computational Molecular and Materials Science, Tel Aviv University, Tel Aviv, Israel 69978

1 Introduction

The Anderson impurity model 1 is a fundamental model for studying strongly correlated open quantum systems that appears in many physical situations, including coupled quantum dots in semiconductor heterostructures 2, 3 and in molecular- and nano-electronics.4, 5 It is one of the simplest models of interacting particles and exhibits complex many-body phenomena, such as the Coulomb blockade,6 pair tunneling,7 and the Kondo effect.8, 9 Such nonequilibrium steady-state effects of interacting open quantum systems, continue to present a grand challenge for theory. As such, the Anderson impurity model is often used as a benchmark for developing approximate methods to study many-body physics of interacting particles, in and out of equilibrium.

In recent years, numerous approximate methods have been developed to study transport through nanoscale interacting systems. Among these are quantum master equation approaches and their generalizations,10, 11, 12, 13, 14, 15, 16 approaches based on the nonequilibrium Green’s function formalism,17, 18, 19, 20, 21, 22, 23, 24 and quasi-classical mapping techniques.25, 26 More recently, numerically exact approaches (namely, methods that allow for a systematic convergence of the results) have been proposed that allow for an assessment of the approximate methods in certain regimes of interactions and temperatures. Most notable are real-time path integral methods based on diagrammatic expansions of the hybridization or onsite interactions,27, 28, 29, 30, 31, 32, 33, 34, 35 renormalization group techniques,36, 37, 38 or many-body wavefunction techniques.39 Benchmarks of the various approximation schemes are often limited to the so-called symmetric Anderson model,30, 25, 40, 41, 42 where the empty and fully occupied states of the impurity are degenerate.

In this work, we show that the steady-state current in the Anderson impurity model at the symmetric point coincides exactly with an equivalent noninteracting model, transition. The results presented in this study are derived from two complementary approximate methods: 1) the quantum master equation (QME) approach and 2) the equation of motion (EOM) nonequilibrium Green’s function (NEGF) approach.43, 44 The QME approach is adequate in the weak system-leads coupling limit and at high temperatures, for arbitrary onsite interactions, whereas the EOM-NEGF approach used here is accurate for small onsite interactions, but is not limited to weak system-bath couplings. These methods do not account for the Kondo effect or pair-tunneling. As such, this work focuses on the many-body physics of the Anderson model, which is manifested in the Coulomb blockade regime.

2 Model and methods

The Anderson impurity model is defined by the Hamiltonian $H=H_{S}+H_{B}+V$, where

[TABLE]

describes the impurity (or dot), referred to simply as the ‘system Hamiltonian’,

[TABLE]

describes the noninteracting fermionic baths (or leads), and

[TABLE]

describes the hybridization between the system and the leads. In the above, $d_{\sigma}^{\dagger}\ \left(d_{\sigma}\right)$ are the creation (annihilation) operators of an electron on the dot with spin $\sigma=\uparrow,\downarrow$ with one-body energy of $\varepsilon_{\sigma}$, $U$ is the onsite Hubbard interaction, $c_{k\sigma}^{\dagger}\ \left(c_{k\sigma}\right)$ are the creation (annihilation) operators of an electron in mode $k$ of the leads with energy $\varepsilon_{k}$, and $t_{k}$ is the hybridization between the dot and mode $k$ in the lead. The coupling to the quasi-continuous leads is modeled by the leads’ spectral function, $\Gamma_{\ell}\left(\omega\right)$, where $t_{k\in\ell}=\sqrt{\Gamma_{\ell}\left(\varepsilon_{k}\right)\Delta\omega/2\pi}$ is the coupling between the dot and the $k$-th mode of the $\ell=\rm L,R$ bath, and $\Delta\omega$ is the discretization of the leads energy spectrum. Throughout, we take Planck’s constant to be 1, and assume that the coupling to the leads is spin-independent. When a wide-band limit is assumed, $\Gamma_{\ell}\left(\omega\right)=\Gamma_{\ell}$ is taken to be independent of $\omega$.

We will focus on two regimes of the model, the non-interacting regime where $U=0$ and $\varepsilon_{\uparrow}=\varepsilon_{\downarrow}=\varepsilon$, and the symmetric point where $\varepsilon_{\uparrow}=\varepsilon_{\downarrow}=-U/2$. The system is held out of equilibrium by an electric bias by setting the leads’ chemical potentials to be $\mu_{\ell}=\varepsilon_{f}+V_{\ell}$, where $\varepsilon_{f}$ is the Fermi energy of the leads and $V_{\ell}$ is the bias on the $\ell$th lead. A symmetric bias is said to be applied when $\varepsilon_{f}=0$. The many-body energy levels in both regimes of the system Hamiltonian, $H_{S}$, are schematically introduced in Fig. 1. An empty dot, $\left|0\right>$ with energy [math], a singly-occupied dot with an electron in either spin, $\left|\sigma\right>=d^{\dagger}_{\sigma}\left|0\right>$ with energy $\varepsilon$, and a doubly occupied (or full) dot $\left|f\right>=d^{\dagger}_{\uparrow}d^{\dagger}_{\downarrow}\left|0\right>$, with energy $2\varepsilon+U$. We may note that in both cases studied all one-body transitions are given by the same energy difference, which is equal to $\varepsilon$. This is a key feature of the symmetric Anderson model, giving rise to identical one-body transition probabilities as the noninteracting case.

We analyze the model using two methods, a quantum master equation approach and a nonequilibrium Green’s function approach. Although NEGF is typically more general and has a wider regime of validity than the QME approach, we present both results, as the perturbative QME approach facilitates a simple derivation. In the supporting material we discuss the validity of the approximations in more detail and compare the two methods in their appropriate regimes.

2.1 Quantum master equation

The quantum master equation description 45, 46 assumes weak-coupling between the dot and the leads, and is derived using second order perturbation theory in the dot-lead coupling strength. The time scale at which correlations in the leads decay should be smaller than the time scale for the leads to induce a significant change in the dot. This implies that the QME is valid when the temperature of the leads is sufficiently high with respect to the coupling strength (see supporting information for details). It is further assumed that there are no initial correlations between the dot and the leads, and that the leads are in local thermal equilibrium. Since the system Hamiltonian is diagonal in the many-body basis and the dot-lead coupling is bi-linear, the populations are coupled directly and not through the coherences. In other words, the equations of motion for the populations (diagonal terms of the reduced density matrix) and the coherences (non-diagonal terms) are decoupled in the many-body basis. Thus, if the initial state of the system is diagonal in the dot-energy eigen-basis or if we are interested in physical quantities that depend on populations alone, then the QME is reduced to a rate equation for the many-body states:

[TABLE]

Here $p_{0}$ is the probability of occupying the empty dot state $\left|0\right>$ , $p_{\uparrow}$ and $p_{\downarrow}$ are the probabilities of occupying the states $\left|\uparrow\right>$ and $\left|\downarrow\right>$ respectively, and $p_{f}$ is the probability of having a fully occupied dot corresponding to the state $\left|f\right>$. The transition matrix is given by

[TABLE]

where $\gamma^{\ell}_{ij}$ is the transition rate from state $j$ to state $i$ induced by the collective $\ell$ bath. The above rates can be calculated using Fermi’s golden-rule and satisfy local detailed balance such that $\gamma^{\ell}_{ij}=\gamma^{\ell}_{ji}\exp[\beta_{\ell}(\varepsilon_{j}-\varepsilon_{i})]$ where $\beta_{\ell}$ is the inverse temperature of the $\ell$-th lead and $\varepsilon_{i}$ is the energy of the eigenstate $i$.

Assuming symmetry between spin up and spin down, $\varepsilon_{\uparrow}=\varepsilon_{\downarrow}\equiv\varepsilon$, in the wide-band approximation, we can explicitly write

[TABLE]

where $f^{\ell}(\varepsilon-\mu_{\ell})=\left(1+\exp[\beta_{\ell}(\varepsilon-\mu_{\ell})]\right)^{-1}$ is the Fermi-distribution with the chemical potential $\mu_{\ell}$. We note that the results presented in the manuscript are not limited to a wide-band approximation and are valid for any even spectral function, $\Gamma_{\ell}(\varepsilon)=\Gamma_{\ell}(-\varepsilon)$.

The current from the $\ell$ bath is calculated according to

[TABLE]

where $N_{i}$ is the number of electrons in state $i$, and $J_{ij}^{\ell}$ is the probability current from state $j$ to state $i$, due to the coupling to the $\ell$ bath, and $e$ is the charge of the electron.

2.2 Equation of motion nonequilibrium Green’s functions

To support our findings we further investigate the Anderson model at the symmetric point using an EOM approach to NEGF.47, 21, 48, 49, 50 We note that the approach taken here is valid to all orders of $\Gamma$ but is exact to order $\Gamma^{2}U$. Thus, in the limit $U\rightarrow 0$ the NEGF approach recovers the exact noninteracting results for any system-bath coupling strength. The approach does not capture the Kondo and pair-tunneling effects but is capable of describing the Coulomb blockade regime.47, 21, 48, 49, 50 We begin by defining the impurity one-body nonequilibrium contour ordered Green’s function,43, 44, 51

[TABLE]

where $\hat{T}_{c}$ is the contour ordering operator and the average $\langle\cdots\rangle$ should be interpreted as the trace over the many particle Hilbert space with the equilibrium density matrix at $t=-\infty$. The equation of motion for the Green’s functions on the Keldysh contour is obtained from the Heisenberg equation for the time propagation of the operators, leading to an equation of motion which introduces new, higher-order Green’s functions, for which equations of motion are also derived. For interacting Hamiltonians, this method produces a hierarchy of inter-dependent equations for Green’s functions of higher orders. A common procedure is to truncate the equations by choosing an appropriate approximation (or closure), after which the equations are solved self-consistently. In this work we use a closure that is exact for the noninteracting limit, where $U=0$, for any value of the system-leads coupling $\Gamma$, and at any temperature, reproducing the exact resonant level model solution. It also provides a good approximation at sufficiently large values of $U$, in the Coulomb blockade regime.47, 22 A detailed derivation of the closure is provided in the Supplementary Information.

We solve the equations of motion on the Keldysh contour and analytically continued to the real time axis using Langreth rules.52 The integral equations for the retarded and lesser NEGFs in steady state are transformed to the frequency domain, resulting in the following algebraic equations:

[TABLE]

where $\bar{\sigma}$ is the opposite spin to $\sigma$,

[TABLE]

and $G\left(\omega\right)=\int_{-\infty}^{\infty}d\tau\,G\left(\tau\right)e^{i\omega\tau}$, with $\tau=t-t^{\prime}$. In the above, the steady state population of spin $\sigma$ is given by

[TABLE]

and we used the following definitions for the noninteracting Green’s functions:

[TABLE]

with $\delta=0$ or $\delta=U$. The self energy due to the coupling to the leads are given by:53

[TABLE]

The above self energy is identical to the resonant level model self energy, where again $f^{\ell}(\omega)$ is the Fermi distribution with inverse temperature $\beta_{\ell}$ and a symmetric spectral function was assumed, $\Gamma_{\ell}\left(\omega\right)=\Gamma_{\ell}\left(-\omega\right)$. The current at steady-state due to the $\ell$-th bath is given by, 53

[TABLE]

3 Results and discussion

In Fig. 2 we plot the magnitude of the current for the noninteracting (left panels) and interacting (middle panels) models, $I_{{\rm L}}^{(0)}$ and $I_{{\rm L}}^{(U)}$, respectively. The results are plotted as a function of the Fermi energy ($\varepsilon_{f}$) and the bias voltage applied between the leads ($V=\mu_{L}-\mu_{R}$), with the chemical potentials of the leads set to $\mu_{{}_{\rm L,R}}=\varepsilon_{f}\pm V/2$. Two different temperatures were considered: $T_{\rm L}=T_{\rm R}=\Gamma$ (upper panels) and $T_{\rm L}=T_{\rm R}=10\Gamma$ (lower panels). To directly compare the results between the interacting and noninteracting models, we also plot the absolute value of the difference between the two (right panels). Results are shown for the EOM-NEGF approach, but a similar qualitative picture emerges within the QME formalism.

Focusing on the results for the noninteracting system (Fig. 2, panels (a1) and (a2)), as expected, we find significant values for the current at a finite bias when $\varepsilon_{f}-\varepsilon<|V|/2$, leading to the well-known diamond-like current characteristics, broadened by the temperature $T$. The picture is a bit more evolved for the interacting case (Fig. 2, panels (b1) and (b2)), where a double diamond-like current characteristic shape is observed for $U=-2\varepsilon$. The two-step value of the current results from the well-known Coulomb blockade, where the bias voltage is not sufficiently large to overcome the onsite repulsion, and only one conductance channel is open at intermediate bias voltage values. Only when $|V|$ becomes sufficiently large compared to $U$ an additional conducting channel opens up, and consequently the current increases to its maximal value. We note that the Coulomb blockade phenomenon is suppressed in the high temperature limit, as would be expected, where the current rises gradually to its maximum value.

Interestingly, when the bias voltage is applied symmetrically, i.e. $\varepsilon_{f}=0$, we find that the current-voltage characteristics are identical for both interacting and noninteracting models, suggesting that the Coulomb blockade is completely suppressed. This is clearly depicted in panels (c1) and (c2) of Fig. 2, where we plot the difference between the interacting and noninteracting currents, which diminishes as we approach the line $\varepsilon_{f}=0$ for all values of $V$. This is clarified in Fig. 3 where we show cuts through the 2D plots of Fig. 2 for two values of $\varepsilon_{f}$. When $\varepsilon_{f}\neq 0$, the I-V curves differ at both intermediate and high temperatures, while for a symmetric bias, $\varepsilon_{f}=0$, the interacting and noninteracting results overlap irrespective of the value of $T$. The results for $\varepsilon_{f}=0$ for different values of $U$ and $\varepsilon$ are summarized in Fig. 4, where we plot the difference between the interacting and noninteracting currents for a fixed bias voltage of $V=4\Gamma$. The green dashed line indicates the symmetric line where $U=-2\varepsilon$, for which the difference between the currents vanishes.

We now turn to rationalize this result. We begin by using the QME approach, which provides a simple relation between the interacting and noninteracting currents at $\varepsilon_{f}=0$. In the case where $\mu_{\rm L}=-\mu_{\rm R}\equiv\mu$ the master equation has a rather simple exact solution for the steady-state populations and currents (see supporting information). The steady-state current is given by

[TABLE]

for the noninteracting case, whereas for the symmetric interacting model, the current takes a slightly more complicated form:

[TABLE]

When the temperatures of both leads are balanced, i.e. $T_{\rm L}=T_{\rm R}$, Eqs. (15) and (16) coincide exactly. Note that the leads’ spectral functions, $\Gamma_{\rm L}$ and $\Gamma_{\rm R}$, need not be equal and the two expressions coincide as long as $\Gamma_{L}$ and $\Gamma_{R}$ are even functions in energy. In the case where $\Gamma_{\rm L}=\Gamma_{\rm R}$, we obtain half-filling of the dot, i.e., n_{\uparrow}\big{|}_{U=-2\varepsilon}=n_{\downarrow}\big{|}_{U=-2\varepsilon}=\frac{1}{2}. Furthermore, the probabilities of the symmetric model are directly related to the non-interacting system average populations by:

[TABLE]

One additional consequence of having the same electric currents is that the steady-state entropy production will also coincide exactly for the two models at this special symmetric point for $\varepsilon_{f}=0$.

The equivalence between the interacting and noninteracting currents can also be derived within the EOM-NEGF approach. For clarity, we restrict the discussion to the wide band limit and assume that $\Gamma_{\rm L}=\Gamma_{\rm R}=\Gamma/2$. In this limit, we find that the difference between the steady state currents of the interacting and noninteracting models is given by (see supporting information):

[TABLE]

When the electric bias is centred around $\varepsilon_{f}=0$ with $T_{\rm L}=T_{\rm R}$, the above difference vanishes exactly, as the integrand becomes an odd function of $\omega$.

4 Conclusions

In conclusion, we have shown that the current in the Coulomb blockade regime for a symmetric Anderson impurity model is identical to the current in a non-interacting model. This was demonstrated (numerically and analytically) within the quantum master equation approach as well as the equation of motion NEGF approach (within a two-particle closure) for a wide range of model parameters, bias voltages, and temperatures. The Anderson impurity model is the canonical model to study weakly and strongly correlated effects away from equilibrium and is routinely used to assess the accuracy of approximate methods to compute the dynamics and steady-state properties. Limiting such assessments to the special symmetric point, as described in this work, is by no means a signature of a valid many-body approximation. Such methods should always be tested away from the symmetric case.

{acknowledgement}

The authors would like to thank Wenjie Dou and Guy Cohen for fruitful discussions. M.T. thanks Eran Rabani for kindly hosting his sabbatical stay at the Chemistry Department of the University of California at Berkeley and acknowledges support by the German Research Foundation (DFG). J.O. acknowledges support from the Georg H. Endress foundation. This work was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Materials Sciences and Engineering Division, under Contract No. DEAC02-05-CH11231 within the Physical Chemistry of Inorganic Nanostructures Program (KC3103).

{suppinfo}

In this supporting information we give a detailed derivation of the master equation and the equation of motion approach to the nonequilibrium Green’s functions for the Anderson impurity model and assess the validity of these approximations.

The Anderson impurity model is defined by the Hamiltonian $H=H_{S}+H_{B}+V$, where

[TABLE]

describes the impurity (or dot), referred to simply as the ‘system Hamiltonian’,

[TABLE]

describes the noninteracting fermionic baths (or leads), and

[TABLE]

describes the hybridization between the system and the leads. In the above, $d_{\sigma}^{\dagger}\ \left(d_{\sigma}\right)$ are the creation (annihilation) operators of an electron on the dot with spin $\sigma=\uparrow,\downarrow$ with one-body energy of $\varepsilon_{\sigma}$, $U$ is the on-site Hubbard interaction, $c_{k\sigma}^{\dagger}\ \left(c_{k\sigma}\right)$ are the creation (annihilation) operators of an electron in mode $k$ of the leads with energy $\varepsilon_{k}$, and $t_{k}$ is the hybridization between the dot and mode $k$ in the lead. The coupling to the quasi-continuous leads is modeled by the leads’ spectral function, $\Gamma_{\ell}\left(\omega\right)$, where $t_{k\in\ell}=\sqrt{\Gamma_{\ell}\left(\varepsilon_{k}\right)\Delta\omega/2\pi}$ is the coupling between the dot and the $k$th mode of the $\ell=L,R$ bath, and $\Delta\omega$ is the discretization of the leads energy spectrum. Throughout, we take Planck’s constant and the Boltzmann factor to be 1.

5 Master equation approach

For weak coupling to the leads and at high temperatures, the master equation (ME) approach is an adequate approximation to describe transport in the Anderson impurity model. In the many-body state basis of the system, the ME reduces to the rate equation:

[TABLE]

where $p_{0}$ is the probability of occupying the empty dot state $\left|0\right>$ , $p_{\uparrow}$ and $p_{\downarrow}$ are the probabilities of occupying a single electron in the states $\left|\uparrow\right>$ and $\left|\downarrow\right>$ respectively, and $p_{f}$ is the probability of having a fully occupied dot corresponding to the state $\left|f\right>$. The transition matrix is given by :

[TABLE]

Here $\gamma^{\ell}_{ij}$ is the transition rate from state $j$ to state $i$ induced by the $\ell=L,R$ bath, which can be calculated using Fermi’s golden rule. Assuming symmetry between spin up and spin down and in the wide-band approximation, the rates are given explicitly by:

[TABLE]

where $\mu_{\ell}$ is the chemical potential for the left ($\ell=L$) or right ($\ell=R$) leads.

The electric, the energy and the heat currents from the $\ell$-bath are given by

[TABLE]

where $\quad j_{ij}^{\ell}=\gamma_{ij}^{\ell}p_{j}(t)-\gamma_{ji}^{\ell}p_{i}(t)$. Here $N_{i}$ and $E_{i}$ are the number of electrons and the energy of the many-body state $i$, and $j_{ij}^{\ell}$ is the probability current from state $j$ to state $i$ induced by the coupling to the $\ell$ bath. The probability currents can then be expressed as:

[TABLE]

Note that if the left and right temperature of the baths are equal $(T_{L}=T_{R})$ and the electric currents of the two systems are the same, then the entropy production $\Delta S$ of the two systems at steady-state is also identical. At steady-state $\Delta S=-\sum_{\ell}\frac{I_{\ell}^{H}}{T_{\ell}}$. Since $T_{L}=T_{R}$ the entropy production depends solely on the electric current.

In the following we assume $\mu_{L}=-\mu_{R}\equiv\mu$.

5.1 Noninteracting model, $U=0$:

Solving equations (22)-(24) for the noninteracting system, the steady-state probabilities are given by:

[TABLE]

Since spin up and down are independent, the single particle probabilities are simply $q_{0}=\sqrt{p_{0}}$ of having no electron and $q_{1}=\sqrt{p_{f}}$ of having a single electron. Thus, the occupation number (average population) is given by,

[TABLE]

The steady-state current can be calculated using equations (25) and (26), and is given by

[TABLE]

5.2 Interacting model, symmetric Anderson model $U=-2\varepsilon$:

The steady state probabilities from equations (22)-(24) are

[TABLE]

The additional necessary condition for having half-filling is that $\Gamma_{L}=\Gamma_{R}$ and $T_{L}=T_{R}$. In this case, Eq. (30) can be simplified to yield:

[TABLE]

and the occupation number is $n_{\uparrow}=n_{\downarrow}=\frac{1}{2}$. The probabilities $q_{0}$ and $q_{1}$ of the noninteracting system are exactly twice the probabilities of the symmetric interacting system, i.e., $q_{0}=2p_{0}=2p_{f}$ and $q_{1}=2p_{\uparrow}=2p_{\downarrow}$.

Using equations (25), (26) and (30) the steady-state current from the left lead is,

[TABLE]

For $T_{L}=T_{R}\equiv T$ the steady-state current takes the form

[TABLE]

which is identical to the current in the noninteracting model Eq. (29).

In summary, the conditions for having equal currents for the interacting and noninteracting systems are: equal temperatures of the left and right leads $T_{L}=T_{R}$, Fermi energy equal to zero $\mu_{L}=-\mu_{R}$ (symmetric distribution of the bias) and $U=-2\varepsilon$. Note that the leads’ spectral functions, $\Gamma_{L}$ and $\Gamma_{R}$, are not necessarily equal as required for obtaining half-filling, but must be an even function in energy.

As a side note we mention that the difference between the currents in the symmetric and noninteracting models may be expanded to first order in the temperature difference between the leads, $\Delta T=T_{R}-T_{L}$, to obtain

[TABLE]

Thus, one may use the noninteracting result to obtain the interacting current in the vicinity of the symmetry point.

6 Nonequilibrium Green’s function approach

We begin by defining the impurity one-body nonequilibrium Green’s function (GF) on the Keldysh contour, 43, 44

[TABLE]

where $\hat{T}_{c}$ is the contour ordering operator. Using the equation of motion (EOM) technique 54, 55, 48, we take the derivative of the GF according to one of the contour times, producing the following equation

[TABLE]

where we defined

[TABLE]

The above satisfy the following equations of motion

[TABLE]

where we also defined

[TABLE]

Since the Anderson model is not analytically solvable, this ever-growing set of inter-dependent equations will not come to a close, and so, we choose an appropriate closure, which is known to be a good approximation in the Coulomb blockade regime. Within this closure, the equation of motion for $F_{k\sigma}^{\left(2\right)}\left(\tau,\tau^{\prime}\right)$ is taken exactly,

[TABLE]

whereas $G_{\sigma}^{\left(3\right)}\left(\tau,\tau^{\prime}\right)$ is treated approximately. If one derives the EOM for $G_{\sigma}^{\left(2\right)}\left(\tau,\tau^{\prime}\right)$ under the noninteracting, uncoupled Hamiltonian, the resulting equation of motion is found to be identical to that of $G_{\sigma}^{\left(3\right)}\left(\tau,\tau^{\prime}\right)$ under the same Hamiltonian (and has the same initial value). Therefore, within this closure, we will use the approximation:

[TABLE]

In this derivation, we will also consider the neglected terms from the closure, in order to keep the full description in mind. Thus, we define

[TABLE]

where $\Gamma F_{\sigma}^{\left({\rm corr}\right)}$ is the difference between the exact $G_{\sigma}^{\left(3\right)}$ and our approximation and the leading contribution due to the coupling to the leads is factored out. As a side note, we can show that the correction term will contribute at least $\propto\Gamma$ by writing out the full equation of motion for $G_{\sigma}^{\left(3\right)}$,

[TABLE]

To arrive at the uncoupled equation, we neglect all terms with $\sum_{k\in L,R}t_{k}$, keeping only terms that appear in the propagation under the dot Hamiltonian,

[TABLE]

thus,

[TABLE]

Solving the equations of motion on the Keldysh contour, one arrives at the following integral equations

[TABLE]

where the noninteracting Green’s function on the contour are defined as

[TABLE]

and the self energy due to the coupling to the leads (identical to the noninteracting model self energy) is given by

[TABLE]

By means of analytical continuation (Langreth rules)52, we find that the equations for the retarded and lesser Green’s functions are given by

[TABLE]

where the noninteracting Green’s functions are given by

[TABLE]

In steady state, all functions depend on time differences rather than two time variables, and thus, the integral equation are simplified to algebraic equations in frequency space. Therefore,

[TABLE]

where the $\sigma$-spin electron steady state population, $\left\langle n_{\sigma}\right\rangle$ is given by:

[TABLE]

and the noninteracting Green’s functions in frequency are given by

[TABLE]

In the above, we assumed that the spectral function characterizing the dot-lead coupling is symmetric, i.e. $\Gamma\left(\omega\right)=\Gamma\left(-\omega\right)$. Note we also set the Fermi energy of the leads to $\varepsilon_{f}=0$, and consider spectral functions that are symmetric around $\varepsilon_{f}$. This will be critical for considering the symmetric point of the Anderson model. The above equations may now be solved analytically.

6.1 Symmetric case - steady state solution

Let us now consider the more special case of the symmetric Anderson model, where:

[TABLE]

The noninteracting GFs in steady state take the form:

[TABLE]

Solving the steady state equations in frequency space, we start with the retarded GF. Rewriting the first line of equation 56, we find

[TABLE]

where we recognized the retarded GF for the noninteracting model 53 with $\varepsilon=-\frac{U}{2}$,

[TABLE]

Similarly, from the second line of 56, we find:

[TABLE]

where

[TABLE]

is the retarded GF for the noninteracting model with $\varepsilon=+\frac{U}{2}$.

Plugging the above into equation (62), we find:

[TABLE]

where we defined the closure-correction contribution to $G^{r}$:

[TABLE]

Finally, using

[TABLE]

we find

[TABLE]

where

[TABLE]

and we took the wide-band limit (WBL) for the lead-dot coupling, $\Sigma^{r}\left(\omega\right)=-\frac{i}{2}\left(\Gamma_{L}+\Gamma_{R}\right)=-\frac{i\Gamma}{2}$.

Next, we look at the lesser GF’s from equation (56):

[TABLE]

Terms of the form $g^{<}\left(1-g^{r}\Sigma^{r}\right)^{-1}$ vanish 22 and we are left with:

[TABLE]

By using $G^{a}=\left(G^{r}\right)^{*}$ and equation (63), we find

[TABLE]

where the lesser GF for the noninteracting model is given by

[TABLE]

and we defined

[TABLE]

After simplifying and taking the wide-band limit for the lesser self energy,

[TABLE]

the difference between the interacting and noninteracting GFs takes the form:

[TABLE]

and the correction term is

[TABLE]

To recap, we have solved the steady state Green’s functions for the symmetric Anderson impurity model, $G_{\sigma}^{r,<}\left(\omega\right)$, where $U=-2\varepsilon_{\uparrow}=-2\varepsilon_{\downarrow}$. These are expressed in terms of the noninteracting model with $\varepsilon=\mp U/2$, described by the Green’s functions $G_{\sigma}^{(0)r,<}\left(\omega\right)$ and $G_{\sigma}^{(0^{\prime})r,<}\left(\omega\right)$, respectively. Correction terms that stem from the closure of the equations of motion, $F^{\left({\rm corr}\right)r,<}\left(\omega\right)$, are also included, to give:

[TABLE]

where

[TABLE]

More explicitly, in the wide-band limit, the equations take the form

[TABLE]

6.2 Steady state current calculation

First recall that in the symmetric case of the Anderson model,

[TABLE]

therefore, we will look at the current of one spin, where the full current will be twice the former. The particle current from the $\ell$ lead in steady state for both the symmetric and the noninteracting models, is given by:

[TABLE]

where the electric current will be given by multiplying by the charge of the electron, $e$. The difference between the Anderson current (for one of the spins) and the noninteracting model current is therefore given by:

[TABLE]

Using the results above, we find:

[TABLE]

where we defined the contribution to the current as a result of the corrections $G_{\sigma}^{\left({\rm corr}\right)}$ due to higher order GFs which were neglected in the closure, as:

[TABLE]

Since $G^{\left({\rm corr}\right)}\sim\Gamma U^{2}$ and $\Gamma^{2}U^{2}$, the correction to the current will contribute as $I^{\left({\rm corr}\right)}\sim\Gamma^{2}U^{2}$ and $\Gamma^{3}U^{2}$.

When the electric bias is centred around the Fermi level, i.e. $\mu_{L}=-\mu_{R}=\mu$, the difference between the noninteracting model and symmetric-Anderson model currents is zero up to the high order correction due to the closure.

[TABLE]

This results shows that within the closure specified above, which captures the Coulomb blockade, at the symmetric point, the noninteracting and interacting currents are equal. Furthermore, the leading corrections to this result scale as $\Gamma^{2}U^{2}$.

7 Comparing the NEGF and ME approaches

To justify the validity of the results obtained using the ME approach we assess the validity of the ME by comparing the results to the NEGF which becomes exact for the noninteracting system. Fig. 5 demonstrates the regime of validity of the ME. We plot the relative error in the steady-state current for the noninteracting system, $|I^{exact}-I^{ME}|/I^{exact}$. The exact current is calculated using the NEGF approach which is exact for $U=0$. We find that when $T,\varepsilon\gtrsim 10\Gamma$ the relative error is $\sim 1\%$.

In Fig. 6 we plot the steady state current for the interacting system using the NEGF and its absolute value difference from the current obtained from the ME, $|I^{NEGF}-I^{ME}|$. At sufficiently high temperatures ($T\sim 10\Gamma$) the currents calculated from NEGF and the ME are in good agreement.