Symmetry breaking and restoration using the equation-of-motion technique for nonequilibrium quantum impurity models

The fundamental object of nonequilibrium many-body theory, is the contour ordered Green function [22, 23, 39],

$$\begin{eqnarray} &&G({\mathbf{r}}_{2},{\tau }_{2},{\mathbf{r}}_{1},{\tau }_{1})=-\frac{\mathrm{i}}{\hbar }\langle {T}_{\mathrm{C}}{\hat {\Psi }}_{\mathrm{H}}({\mathbf{r}}_{2},{\tau }_{2}){\hat {\Psi }}_{\mathrm{H}}^{\dagger }({\mathbf{r}}_{1},{\tau }_{1})\rangle , \end{eqnarray} \tag{ 3 }$$

where the time variables (τ₁ and τ₂) are along the Schwinger–Keldysh contour as depicted in figure 1, T_C is the contour time ordering operator which moves operators with earlier time arguments to the right, ${\hat {\Psi }}_{\mathrm{H}}/{\hat {\Psi }}_{\mathrm{H}}^{\dagger }$ are the system's annihilation/creation field operators in the Heisenberg picture (in what follows we omit the H index), and the average 〈⋯〉 should be interpreted as the trace over the many particle Hilbert space with the equilibrium (t =− ∞) density matrix. The EOM for the contour ordered GF is obtained from the Heisenberg EOM for a Heisenberg operator $\frac{\mathrm{d}}{\mathrm{d}t}\hat {A}(t)=\frac{\mathrm{i}}{\hbar }[\hat {H}(t),\hat {A}(t)]+\frac{\partial }{\partial t}\hat {A}(t)$, where in our case $\hat {H}(t)={\hat {H}}_{0}+\hat {V}(t)$. Here ${\hat {H}}_{0}$ stands for the one-body non-interacting part of $\hat {H}(t),\hat {V}(t)=\hat {H}(t)-{\hat {H}}_{0}$, and $[\hat {A},\hat {B}]$ is the commutator. Let us consider a generic example, The EOM [40] for G(r₂,τ₂,r₁,τ₁) can be written as (omitting the r dependence for brevity)

$$\begin{eqnarray} G({\tau }_{2},{\tau }_{1})&=&{g}_{2}({\tau }_{2},{\tau }_{1})\langle \{ \hat {\Psi },{\hat {\Psi }}^{\dagger }\} \rangle -\frac{\mathrm{i}}{\hbar }\int _{\mathrm{C}}\mathrm{d}\tau \hspace{0.167em} {g}_{2}({\tau }_{2},\tau )\nonumber\\ &&\times \, \langle {T}_{\mathrm{C}}[\hat {\Psi }(\tau ),\hat {V}(\tau )]{\hat {\Psi }}^{\dagger }({\tau }_{1})\rangle , \end{eqnarray} \tag{ 4 }$$

where $(\mathrm{i}\hbar \frac{\partial }{\partial {\tau }_{2}}-\varepsilon ){g}_{2}({\tau }_{2},{\tau }_{1})=\delta ({\tau }_{1}-{\tau }_{2}),\{ \hat {A},\hat {B}\}$ is the anti-commutator, and ε is defined from the equation, $\varepsilon \hat {\Psi }(\tau )=[\hat {\Psi }(\tau ),{\hat {H}}_{0}]$. For example, if ${\hat {H}}_{0}={\sum }_{n}({\varepsilon }_{n}-\mu ){\hat {d}}_{n}^{\dagger }{\hat {d}}_{n}$ and $\hat {\Psi }={\hat {d}}_{i}$ then ε = ε_i − μ. Following Langreth theorem [41], we can change the contour integration in equation (4) to integration along the real-time axis. This yields (see section 3 for the definitions of the different real-time GFs)

$$\begin{eqnarray} \fl {G}^{r}({t}_{2},{t}_{1})&=&{g}_{2}^{r}({t}_{2},{t}_{1})\langle \{ \hat {\Psi },{\hat {\Psi }}^{\dagger }\} \rangle \nonumber\\ &&+~\int \nolimits _{{t}_{1}}^{{t}_{2}}\mathrm{d}t \hspace{0.167em} {g}_{2}^{r}({t}_{2},t){\mathbb{G}}^{r}(t,{t}_{1}), \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \fl {G}^{\lt }({t}_{2},{t}_{1})&=&{g}_{2}^{\lt }({t}_{2},{t}_{1})\langle \{ \hat {\Psi },{\hat {\Psi }}^{\dagger }\} \rangle \nonumber\\ &&+\int \nolimits _{{t}_{0}}^{{t}_{2}}\mathrm{d}t \hspace{0.167em} {g}_{2}^{r}({t}_{2},t){\mathbb{G}}^{\lt }(t,{t}_{1})\nonumber\\ &&+\int \nolimits _{{t}_{0}}^{{t}_{1}}\mathrm{d}t \hspace{0.167em} {g}_{2}^{\lt }({t}_{2},t){\mathbb{G}}^{a}(t,{t}_{1}), \end{eqnarray} \tag{ 6 }$$

where the time variables (t,t₁ and t₂) are along the real-time axis. G^r(t₂,t₁) is the retarded GF usually used to calculate the response of the system at time t₂ to an earlier perturbation of the system at time t₁. G^<(t₂,t₁) is the lesser GF which plays the role of the single particle density matrix, and $\mathbb{G}({\tau }_{2},{\tau }_{1})=-\frac{\mathrm{i}}{\hbar }\langle {T}_{\mathrm{C}}[\hat {\Psi }({\tau }_{2}),\hat {V}({\tau }_{2})]{\hat {\Psi }}^{\dagger }({\tau }_{1})\rangle$ is a new GF generated by the EOM procedure. Depending on the Hamiltonian it can be a single particle GF or a many particle GF and can involve lead operators as well as system operators. In steady state, the GFs depend only on the difference in time, t = t₂ − t₁, which is simpler to express in Fourier space

$$\begin{eqnarray} {G}^{r}(\omega )&=&{g}_{2}^{r}(\omega )\langle \{ \hat {\Psi },{\hat {\Psi }}^{\dagger }\} \rangle +{g}_{2}^{r}(\omega ){\mathbb{G}}^{r}(\omega ), \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} {G}^{\lt }(\omega )&=&{g}_{2}^{\lt }(\omega )\langle \{ \hat {\Psi },{\hat {\Psi }}^{\dagger }\} \rangle +{g}_{2}^{r}(\omega ){\mathbb{G}}^{\lt }(\omega )\nonumber\\ &&+~{g}_{2}^{\lt }(\omega ){\mathbb{G}}^{a}(\omega ). \end{eqnarray} \tag{ 8 }$$

To simplify the notation we denote the Fourier transform of G(t₂ − t₁) = G(t) as G(ω), i.e., functions with an argument 'ω' are Fourier transforms of their time-domain counterparts. At this stage one has to evaluate G(t₂,t₁) (G(t) in steady state). Except for very simple cases, where an exact closure can be obtained, writing the EOM for G(t₂,t₁) will produce new and/or 'higher order' GFs that need to be evaluated. This leads (in principle) to an infinite set of equations. The idea of the EOM method is therefore, to truncate this hierarchy of equations making a mean-field like approximation for the 'higher order' GFs through lower order functions. This is the Achilles heel of this method as there is no systematic way to close the equations. Usually the approximations have physical meaning within the regime of the problem at hand [34, 42, 43]. In what follows we demonstrate that different approximations can sometimes break symmetry relations that the GFs must fulfil. We will use two impurity models to demonstrate at what level of approximation the symmetry relations are violated and propose a scheme to restore symmetrization.

Zoom In Zoom Out Reset image size

To illustrate the shortcomings of the EOM approach, we refer to the Anderson model [34, 37, 44] and the double Anderson model [38] to represent two different degrees of complexity in correlated systems. As commonly used, we split the total Hamiltonian into three parts [28]:

$$\begin{eqnarray} &&\hat {H}={\hat {H}}_{\mathrm{sys}}+{\hat {H}}_{\mathrm{bath}}+{\hat {H}}_{\mathrm{int}}, \end{eqnarray} \tag{ 9 }$$

where ${\hat {H}}_{\mathrm{bath}}$ describes the macroscopic leads (left and right contacts), ${\hat {H}}_{\mathrm{sys}}$ describes the system of interest (in our case the impurities), and ${\hat {H}}_{\mathrm{int}}$ is the interaction Hamiltonian between the system and the leads. The contacts (leads) are modelled as infinite non-interacting fermionic baths [45–47] with a Hamiltonian in second quantization given by

$$\begin{eqnarray} &&{\hat {H}}_{\mathrm{bath}}=\mathop{\sum }\limits _{\sigma ,k\in \{ \mathrm{L},\mathrm{R}\} }{\epsilon }_{k \sigma }c_{k \sigma }^{\dagger }{c}_{k \sigma }, \end{eqnarray} \tag{ 10 }$$

where _kσ is the energy of a free electron in the left (L) or right (R) lead, in momentum state k and spin σ. The operator ${c}_{k \sigma }\hspace{0.167em} ({c}_{k \sigma }^{\dagger })$ is the annihilation (creation) operator of such an electron. The form chosen for ${\hat {H}}_{\mathrm{sys}}$ depends on the system studied. For the Anderson impurity model [37]

$$\begin{eqnarray} &&{\hat {H}}_{\mathrm{sys}}=\mathop{\sum }\limits _{\sigma \in \{ \uparrow ,\downarrow \} }{\epsilon }_{\sigma }{n}_{\sigma }+U{n}_{\uparrow }{n}_{\downarrow }. \end{eqnarray} \tag{ 11 }$$

Here ${n}_{\sigma }={d}_{\sigma }^{\dagger }{d}_{\sigma }$ is the number operator of the spin σ electron with energy ε_σ and U is the repulsion energy between two electrons on the same site with opposite spins (intra-site repulsion). The second model we discuss is the double Anderson model [38]

$$\begin{eqnarray} {\hat {H}}_{\mathrm{sys}}&=&\mathop{\sum }\limits _{\sigma ,m\in \{ \alpha ,\beta \} }{\epsilon }_{m \sigma }{n}_{m \sigma }+\mathop{\sum }\limits _{m}{U}_{m}{n}_{m \uparrow }{n}_{m \downarrow }\nonumber\\ &&+~\mathop{\sum }\limits _{\sigma ,{\sigma }^{\prime}}V_{\alpha \beta }^{\sigma {\sigma }^{\prime}}{n}_{\alpha \sigma }{n}_{\beta {\sigma }^{\prime}}+\mathop{\sum }\limits _{\sigma }[h_{\alpha \beta }^{\sigma }d_{\alpha }^{\dagger }{d}_{\beta }+\text{h.c}.],\nonumber\\ \end{eqnarray} \tag{ 12 }$$

where the first two terms on the RHS are similar to the Anderson impurity model Hamiltonian (extended to 2 sites), ${V}_{\alpha \beta }^{\sigma {\sigma }^{\prime}}$ is the repulsion energy between two electrons on different sites (inter-site repulsion), and ${h}_{\alpha \beta }^{\sigma }$ is the coupling strength for electron hopping between the two sites. The interaction between the system and the contacts is simply given by the tunnelling Hamiltonian [48]

$$\begin{eqnarray} &&{\hat {H}}_{\mathrm{int}}=\mathop{\sum }\limits _{m,\sigma ,k\in \{ \mathrm{L},\mathrm{R}\} }t_{k m}^{\sigma }{c}_{k \sigma }^{\dagger }{d}_{m \sigma }+\text{h.c.} \end{eqnarray} \tag{ 13 }$$

The parameter ${t}_{k m}^{\sigma }$ represents the coupling strength between the system and the leads, and the index m runs over the site index {α,β} in the double Anderson model.

In the Keldysh formalism the two time NEGF is defined on a contour (figure 1) [22]. In accordance with where on the contour the two times are placed one can define four real-time GFs [49]; the time ordered G^t, anti-time ordered ${G}^{\bar {t}}$, lesser G^<, and greater G^> GF:

$$\begin{eqnarray} &&\fl {G}_{\alpha \beta }({\tau }_{2},{\tau }_{1})\nonumber\\ &&=~\left \{ \begin{array}{@{}l@{}} \displaystyle {G}_{\alpha \beta }^{t}({t}_{2},{t}_{1})=-\frac{\mathrm{i}}{\hbar }\langle T{\Psi }_{\alpha }({t}_{2}){\Psi }_{\beta }^{\dagger }({t}_{1})\rangle \\\tqs \displaystyle {\tau }_{1},{\tau }_{2}\in {C}_{+} \\ \displaystyle {G}_{\alpha \beta }^{\bar {t}}({t}_{2},{t}_{1})=-\frac{\mathrm{i}}{\hbar }\langle \bar {T}{\Psi }_{\alpha }({t}_{2}){\Psi }_{\beta }^{\dagger }({t}_{1})\rangle \\\tqs \displaystyle {\tau }_{1},{\tau }_{2}\in {C}_{-} \\ \displaystyle {G}_{\alpha \beta }^{\lt }({t}_{2},{t}_{1})=+\frac{\mathrm{i}}{\hbar }\langle {\Psi }_{\beta }^{\dagger }({t}_{1}){\Psi }_{\alpha }({t}_{2})\rangle \\\tqs \displaystyle {\tau }_{2}\in {C}_{+},\quad {\tau }_{1}\in {C}_{-} \\ \displaystyle {G}_{\alpha \beta }^{\gt }({t}_{2},{t}_{1})=-\frac{\mathrm{i}}{\hbar }\langle {\Psi }_{\alpha }({t}_{2}){\Psi }_{\beta }^{\dagger }({t}_{1})\rangle \\\tqs \displaystyle {\tau }_{1}\in {C}_{+},\quad {\tau }_{2}\in {C}_{-}. \end{array}\right . \end{eqnarray} \tag{ 14 }$$

Here, T is the time ordering operator, and $\bar {T}$ is the anti-time ordering operator which orders operators in the opposite way of T. Using the definitions in equation (14), the retarded and advanced GFs are defined as:

$$\begin{eqnarray} {G}_{\alpha \beta }^{r}({t}_{2},{t}_{1})&=&{G}_{\alpha \beta }^{t}({t}_{2},{t}_{1})-{G}_{\alpha \beta }^{\lt }({t}_{2},{t}_{1})\nonumber\\ &=&-\frac{\mathrm{i}}{\hbar }\theta ({t}_{2}-{t}_{1})\langle \{ {\Psi }_{\alpha }({t}_{2}),{\Psi }_{\beta }^{\dagger }({t}_{1})\} \rangle , \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} {G}_{\alpha \beta }^{a}({t}_{2},{t}_{1})&=&{G}_{\alpha \beta }^{\lt }({t}_{2},{t}_{1})-{G}_{\alpha \beta }^{\bar {t}}({t}_{2},{t}_{1})\nonumber\\ &=&\frac{\mathrm{i}}{\hbar }\theta ({t}_{1}-{t}_{2})\langle \{ {\Psi }_{\alpha }({t}_{2}),{\Psi }_{\beta }^{\dagger }({t}_{1})\} \rangle . \end{eqnarray} \tag{ 16 }$$

As already mentioned, the retarded GF can be used to calculate the response of the system at time t₂ to an earlier perturbation of the system at time t₁ and is proportional to the local density of states, while the lesser GF is also known as the particle propagator and plays the role of the single particle density matrix. Using equations (14)–(16) it is clear that the following relations must hold:

$$\begin{eqnarray} \begin{array}{@{}rcl@{}} \displaystyle {G}_{\alpha \beta }^{r}({t}_{2},{t}_{1})&\displaystyle =&\displaystyle ({G}_{\beta \alpha }^{a}({t}_{1},{t}_{2}))^{\ast },\\ \displaystyle {G}_{\alpha \beta }^{\lt ,\gt }({t}_{2},{t}_{1})&\displaystyle =&\displaystyle -({G}_{\beta \alpha }^{\lt ,\gt }({t}_{1},{t}_{2}))^{\ast },\\ \displaystyle {G}_{\alpha \beta }^{r}({t}_{2},{t}_{1})-{G}_{\alpha \beta }^{a}({t}_{2},{t}_{1})&\displaystyle =&\displaystyle {G}_{\alpha \beta }^{\gt }({t}_{2},{t}_{1})-{G}_{\alpha \beta }^{\lt }({t}_{2},{t}_{1}). \end{array} \end{eqnarray} \tag{ 17 }$$

In steady state these relations can be rewritten in Fourier space as:

$$\begin{eqnarray} \begin{array}{@{}rcl@{}} \displaystyle {G}_{\alpha \beta }^{r}(\omega )&\displaystyle =&\displaystyle ({G}_{\beta \alpha }^{a}(\omega ))^{\ast },\\ \displaystyle {G}_{\alpha \beta }^{\lt ,\gt }(\omega )&\displaystyle =&\displaystyle -({G}_{\beta \alpha }^{\lt ,\gt }(\omega ))^{\ast },\\ \displaystyle {G}_{\alpha \beta }^{r}(\omega )-{G}_{\alpha \beta }^{a}(\omega )&\displaystyle =&\displaystyle {G}_{\alpha \beta }^{\gt }(\omega )-{G}_{\alpha \beta }^{\lt }(\omega ). \end{array} \end{eqnarray} \tag{ 18 }$$

In what follows we show that the relations in equation (18) do not always hold when the GFs are obtained by the EOM technique with an arbitrary closure, thus leading to unphysical behaviour [50].

Following the derivation in [28, 34, 48] we define the following contour ordered GF:

$$\begin{eqnarray} {G}_{\sigma \sigma }(\tau ,{\tau }^{\prime})&=&-\frac{\mathrm{i}}{\hbar }\langle {T}_{\mathrm{C}}{d}_{\sigma }(\tau ){d}_{\sigma }^{\dagger }({\tau }^{\prime})\rangle , \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} {G}_{2}(\tau ,{\tau }^{\prime})&=&-\frac{\mathrm{i}}{\hbar }\langle {T}_{\mathrm{C}}{n}_{\bar {\sigma }}(\tau ){d}_{\sigma }(\tau ){d}_{\sigma }^{\dagger }({\tau }^{\prime})\rangle , \end{eqnarray} \tag{ 20 }$$

where $\bar {\sigma }$ is the opposite spin of σ. Various approximate decoupling procedures can be applied to the many particle GF [51]. Here we follow the approximation scheme used in [28, 34] where all electronic correlations containing at most one lead operator, are not decoupled and their EOM are calculated. Higher order GFs involving (opposite) spin correlations in the leads are set to zero, and the remaining higher order GFs involving lead and system degrees of freedom are decoupled such that ${F}_{2}(\tau ,{\tau }^{\prime})=-\frac{\mathrm{i}}{\hbar }\langle {T}_{\mathrm{C}}{c}_{k\bar {\sigma }}^{\dagger }(\tau ){d}_{\sigma }(\tau ){c}_{q\bar {\sigma }}(\tau ){d}_{\sigma }^{\dagger }({\tau }^{\prime})\rangle =-{\delta }_{k q}{f}_{\mathrm{L}/\mathrm{R}}({\varepsilon }_{k\bar {\sigma }}-{\mu }_{\mathrm{L}/\mathrm{R}}){G}_{\sigma \sigma }(\tau ,{\tau }^{\prime})$. Given this set of approximations the EOMs for G_σσ(τ,τ') can now close.

In what follows, we show that the symmetry relation for the lesser GF, ${G}_{\sigma \sigma }^{\lt }(\omega )=-({G}_{\sigma \sigma }^{\lt }(\omega ))^{\ast }$, is not satisfied. In turn, this implies that $\langle {n}_{\sigma }\rangle =-\frac{\mathrm{i}\hbar }{2 \pi }\int \nolimits _{-\infty }^{\infty }{G}_{\sigma \sigma }^{\lt }(\omega )\hspace{0.167em} \mathrm{d}\omega$ (the occupation number) is a complex number, which of course is not physical. Following the same derivation given below one can show that the greater GF is not an imaginary function either. All other symmetry requirements are fulfilled for the Anderson impurity model.

Within the approximation scheme of [34], the lesser GF is given by:

$$\begin{eqnarray} &&\fl {G}_{\sigma \sigma }^{\lt }(\omega )\nonumber\\ &&=~{g}^{\lt }(\omega )-{g}^{\lt }(\omega )U{P}^{a}(\omega ){g}_{2}^{a}(\omega ){\Sigma }_{1}^{a}(\omega ){g}^{a}(\omega )\nonumber\\ &&-~{g}^{r}(\omega )U{P}^{r}(\omega ){g}_{2}^{\lt }(\omega ){\Sigma }_{1}^{a}(\omega ){g}^{a}(\omega )\nonumber\\ &&-~{g}^{r}(\omega )U{P}^{r}(\omega ){g}_{2}^{r}(\omega )({\Sigma }_{1}^{r}(\omega ){g}^{\lt }(\omega )+{\Sigma }_{1}^{\lt }(\omega ){g}^{a}(\omega ))\nonumber\\ &&-~{g}^{r}(\omega )U{P}^{r}(\omega ){g}_{2}^{r}(\omega ){\Sigma }_{1}^{r}(\omega ){g}^{\lt }(\omega )\nonumber\\ &&\times ~U{P}^{a}(\omega ){g}_{2}^{a}(\omega ){\Sigma }_{1}^{a}(\omega ){g}^{a}(\omega )\nonumber\\ &&-~{g}^{r}(\omega )U{P}^{r}(\omega ){g}_{2}^{r}(\omega ){\Sigma }_{1}^{\lt }(\omega ){g}^{a}(\omega )\nonumber\\ &&\times ~U{P}^{a}(\omega ){g}_{2}^{a}(\omega ){\Sigma }_{1}^{a}(\omega ){g}^{a}(\omega )\nonumber\\ &&-~{g}^{r}(\omega )U{P}^{r}(\omega ){g}_{2}^{\lt }(\omega ){\Sigma }_{1}^{a}(\omega ){g}^{a}(\omega )\nonumber\\ &&\times ~U{P}^{a}(\omega ){g}_{2}^{a}(\omega ){\Sigma }_{1}^{a}(\omega ){g}^{a}(\omega )\nonumber\\ &&+~U\langle {n}_{\bar {\sigma }}\rangle [{g}^{r}(\omega ){P}^{r}(\omega ){g}_{2}^{\lt }(\omega )+{g}^{\lt }(\omega ){P}^{a}(\omega ){g}_{2}^{a}(\omega )\nonumber\\ &&-~{g}^{r}(\omega ){P}^{r}(\omega ){g}_{2}^{r}(\omega ){\Sigma }_{1}^{r}(\omega ){g}^{\lt }(\omega )U{P}^{a}(\omega ){g}_{2}^{a}(\omega )\nonumber\\ &&-~{g}^{r}(\omega ){P}^{r}(\omega ){g}_{2}^{r}(\omega ){\Sigma }_{1}^{\lt }(\omega ){g}^{a}(\omega )U{P}^{a}(\omega ){g}_{2}^{a}(\omega )\nonumber\\ &&-~{g}^{r}(\omega ){P}^{r}(\omega ){g}_{2}^{\lt }(\omega ){\Sigma }_{1}^{a}(\omega ){g}^{a}(\omega )U{P}^{a}(\omega ){g}_{2}^{a}(\omega )], \end{eqnarray} \tag{ 21 }$$

with

$$\begin{eqnarray} {g}^{r,a}(\omega )&=&(\hbar \omega \pm \mathrm{i}{0}^{+}-{\varepsilon }_{\sigma }-{\Sigma }_{0}^{r,a}(\omega ))^{-1}, \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} {g}_{2}^{r,a}(\omega )&=&(\hbar \omega \pm \mathrm{i}{0}^{+}-{\varepsilon }_{\sigma }-U-{\Sigma }_{4}^{r,a}(\omega ))^{-1}, \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} {g}^{\lt }(\omega )&=&{g}^{r}(\omega ){\Sigma }_{0}^{\lt }(\omega ){g}^{a}(\omega ), \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} {g}_{2}^{\lt }(\omega )&=&{g}_{2}^{r}(\omega ){\Sigma }_{4}^{\lt }(\omega ){g}_{2}^{a}(\omega ), \end{eqnarray} \tag{ 25 }$$

and

$$\begin{eqnarray} &&{P}^{r,a}(\omega )=(1+{g}_{2}^{r,a}(\omega ){\Sigma }_{1}^{r,a}(\omega ){g}^{r,a}(\omega )U)^{-1}. \end{eqnarray} \tag{ 26 }$$

The self-energies are given by:

$$\begin{eqnarray} &&\fl {\Sigma }_{0}^{r,a}(\omega )=\mathop{\sum }\limits _{k\in \{ \mathrm{L,R}\} }\vert {t}_{k \sigma }\vert ^{2}(\hbar \omega -{\varepsilon }_{k \sigma }\pm \mathrm{i}{0}^{+})^{-1}, \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} &&\fl {\Sigma }_{j}^{r,a}(\omega )=\mathop{\sum }\limits _{k\in \{ \mathrm{L,R}\} }A_{k}^{(j)}\vert {t}_{k \sigma }\vert ^{2}[(\hbar \omega +{\varepsilon }_{k\bar {\sigma }}-{\varepsilon }_{\sigma }-{\varepsilon }_{\bar {\sigma }}\nonumber\\ &&-~U\pm \mathrm{i}{0}^{+})^{-1}+(\hbar \omega -{\varepsilon }_{k\bar {\sigma }}-{\varepsilon }_{\sigma }+{\varepsilon }_{\bar {\sigma }}\pm \mathrm{i}{0}^{+})^{-1}],\nonumber\\ &&~j=1,3 \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} &&\fl {\Sigma }_{4}^{r,a,\lt }(\omega )={\Sigma }_{0}^{r,a,\lt }(\omega )+{\Sigma }_{3}^{r,a,\lt }(\omega ). \end{eqnarray} \tag{ 29 }$$

Here ${\Sigma }_{0}^{r,a,\lt }(\omega )$ is the exact self-energy for the non-interacting case, while ${\Sigma }_{1}^{r,a,\lt }(\omega )$ and ${\Sigma }_{3}^{r,a,\lt }(\omega )$ are the self-energies due to the tunnelling of the $\bar {\sigma }$ electron, with ${A}_{k}^{(1)}={f}_{k}({\varepsilon }_{k \sigma }-{\mu }_{k})$ and ${A}_{k}^{(3)}=1$.

The lesser self-energies are defined as in [35]:

$$\begin{eqnarray} {\Sigma }_{j}^{\lt }(\omega )&=&{\Sigma }_{j L}^{\lt }(\omega )+{\Sigma }_{j R}^{\lt }(\omega )\nonumber\\ &=&\mathrm{i}({f}_{\mathrm{ L}}(\varepsilon -{\mu }_{\mathrm{L}}){\Gamma }_{j L}(\omega )\nonumber\\ &&+~{f}_{\mathrm{R}}(\varepsilon -{\mu }_{\mathrm{R}}){\Gamma }_{j R}(\omega )), \end{eqnarray} \tag{ 30 }$$

where ${\Gamma }_{j K}(\omega )=-2\mathrm{Im}({\Sigma }_{j K}^{r}(\omega ))$, j = 0,1,3, and K = L,R.

Applying the principle of reductio ad absurdum, we assume ${G}_{\sigma \sigma }^{\lt }(\omega )$ to be imaginary. Since it must hold for any real value of $\langle {n}_{\bar {\sigma }}\rangle$ between 0 and 1, we argue that the terms multiplying $U\langle {n}_{\bar {\sigma }}\rangle$ in equation (21) must be imaginary. That is (omitting (ω) for brevity)

$$\begin{eqnarray} {A}_{1}&=&{g}^{r}{P}^{r}{g}_{2}^{\lt }+{g}^{\lt }{P}^{a}{g}_{2}^{a}-U[{g}^{r}{P}^{r}{g}_{2}^{r}{\Sigma }_{1}^{r}{g}^{\lt }{P}^{a}{g}_{2}^{a}\nonumber\\ &&+~{g}^{r}{P}^{r}{g}_{2}^{r}{\Sigma }_{1}^{\lt }{g}^{a}U{P}^{a}{g}_{2}^{a}+{g}^{r}{P}^{r}{g}_{2}^{\lt }{\Sigma }_{1}^{a}{g}^{a}U{P}^{a}{g}_{2}^{a}],\nonumber\\ \end{eqnarray} \tag{ 31 }$$

must be imaginary. Moreover, Since A₁ must be imaginary for any value of U, the terms not multiplying it, should be imaginary as well, i.e., 

$$\begin{eqnarray} &&{A}_{2}={g}^{r}{P}^{r}{g}_{2}^{\lt }+{g}^{\lt }{P}^{a}{g}_{2}^{a}, \end{eqnarray} \tag{ 32 }$$

should be imaginary. Up to here we only used the fact that $\langle {n}_{\bar {\sigma }}\rangle$ is a real quantity and we assumed U to be real as well. By definition ${g}_{2}^{\lt }$ and g^< are imaginary, thus, for A₂ to be imaginary, the following identity must hold:

$$\begin{eqnarray} &&\mathrm{Im}({g}^{r}{P}^{r}){g}_{2}^{\lt }+\mathrm{Im}({P}^{a}{g}_{2}^{a}){g}^{\lt }=0. \end{eqnarray} \tag{ 33 }$$

In other words we demand Re(A₂) = 0. Substituting equations (24) and (25) into (33), we get:

$$\begin{eqnarray} &&\mathrm{Im}({g}^{r}{P}^{r}){g}_{2}^{r}{\Sigma }_{4}^{\lt }{g}_{2}^{a}=-\mathrm{Im}({P}^{a}{g}_{2}^{a}){g}^{r}{\Sigma }_{0}^{\lt }{g}^{a}. \end{eqnarray} \tag{ 34 }$$

Thus, if the equality in equation (34) is true, our assumption that ${G}_{\sigma \sigma }^{\lt }$ is imaginary is correct. Using the definitions for g^r,a and ${g}_{2}^{r,a}$ the last equality can be rewritten as:

$$\begin{eqnarray} &&\fl \mathrm{Im}({g}^{r}{P}^{r})\frac{-\mathrm{i}f(\omega )\mathrm{Im}({\Sigma }_{4}^{r})}{(\hbar \omega -{\varepsilon }_{4}-U)^{2}+(\mathrm{Im}({\Sigma }_{4}^{r}))^{2}}\nonumber\\ &&=~\mathrm{Im}({P}^{a}{g}_{2}^{a})\frac{\mathrm{i}f(\omega )\mathrm{Im}({\Sigma }_{0}^{r})}{(\hbar \omega -{\varepsilon }_{0})^{2}+(\mathrm{Im}({\Sigma }_{0}^{r}))^{2}}, \end{eqnarray} \tag{ 35 }$$

where ${\varepsilon }_{0}={\varepsilon }_{\sigma }+\mathrm{Re}({\Sigma }_{0}^{r})$ and ${\varepsilon }_{4}={\varepsilon }_{\sigma }+\mathrm{Re}({\Sigma }_{4}^{r})$. Starting with the LHS of equation (35), we look at g^rP^r:

$$\begin{eqnarray} {g}^{r}{P}^{r}&=&(\hbar \omega -{\varepsilon }_{4}-U-\mathrm{i}\mathrm{Im}({\Sigma }_{4}^{r}))[(\hbar \omega -{\varepsilon }_{0}-\mathrm{i}\mathrm{Im}({\Sigma }_{0}^{r}))\nonumber\\ &&\times \, (\hbar \omega -{\varepsilon }_{4}-U-\mathrm{i}\mathrm{Im}({\Sigma }_{4}^{r}))+U\mathrm{Re}({\Sigma }_{1}^{r})\nonumber\\ &&+~\mathrm{i}U\mathrm{Im}({\Sigma }_{1}^{r})]^{-1}. \end{eqnarray} \tag{ 36 }$$

To simplify the notation we denote: a₀ = ħω − ε₀, ${a}_{1}=\mathrm{Re}({\Sigma }_{1}^{r})$, a₄ = ħω − ε₄ − U and ${b}_{j}=\mathrm{Im}({\Sigma }_{j}^{r})$ with j = 0,1,4. Rewriting the equation for g^rP^r:

$$\begin{eqnarray} &&{g}^{r}{P}^{r}=\frac{{a}_{4}-\mathrm{i}{b}_{4}}{{a}_{0}{a}_{4}-{b}_{0}{b}_{4}+U{a}_{1}-\mathrm{i}({a}_{0}{b}_{4}+{a}_{4}{b}_{0}-U{b}_{1})}.\nonumber\\ \end{eqnarray} \tag{ 37 }$$

Denote D = a₀a₄ − b₀b₄ + Ua₁ and E = a₀b₄ + a₄b₀ − Ub₁, so we can rewrite equation (37) as:

$$\begin{eqnarray} &&{g}^{r}{P}^{r}=\frac{{a}_{4}D+{b}_{4}E-\mathrm{i}({b}_{4}D-{a}_{4}E)}{{D}^{2}+{E}^{2}}. \end{eqnarray} \tag{ 38 }$$

Finally

$$\begin{eqnarray} &&\mathrm{Im}({g}^{r}{P}^{r})=-\frac{{b}_{4}D-{a}_{4}E}{{D}^{2}+{E}^{2}}, \end{eqnarray} \tag{ 39 }$$

and the LHS of equation (35) is given by:

$$\begin{eqnarray} &&\text{LHS}=\frac{\mathrm{i}f(\omega ){b}_{4}}{({a}_{4})^{2}+({b}_{4})^{2}}\frac{{b}_{4}D-{a}_{4}E}{{D}^{2}+{E}^{2}}. \end{eqnarray} \tag{ 40 }$$

Now we turn to analyse the RHS of equation (35). We start by evaluating ${P}^{a}{g}_{2}^{a}$. Using the same notation as we did for g^rP^r, we get:

$$\begin{eqnarray} &&{P}^{a}{g}_{2}^{a}=\frac{{g}_{2}^{a}}{1+{g}_{2}^{a}{\Sigma }_{1}^{a}{g}^{a}U}=\frac{{a}_{0}D+{b}_{0}E+\mathrm{i}({b}_{0}D-{a}_{0}E)}{{D}^{2}+{E}^{2}}.\nonumber\\ \end{eqnarray} \tag{ 41 }$$

Here we used Im(Σ^r) =− Im(Σ^a). Finally:

$$\begin{eqnarray} &&\mathrm{Im}({P}^{a}{g}_{2}^{a})=\frac{{b}_{0}D-{a}_{0}E}{{D}^{2}+{E}^{2}}. \end{eqnarray} \tag{ 42 }$$

The RHS of equation (35) is given by:

$$\begin{eqnarray} &&\text{RHS}=\frac{{b}_{0}D-{a}_{0}E}{{D}^{2}+{E}^{2}}\frac{\mathrm{i}f(\omega ){b}_{0}}{({a}_{0})^{2}+({b}_{0})^{2}}. \end{eqnarray} \tag{ 43 }$$

The equality (equation (35)) now reads:

$$\begin{eqnarray} &&\frac{{b}_{4}}{({a}_{4})^{2}+({b}_{4})^{2}}\frac{{b}_{4}D-{a}_{4}E}{{D}^{2}+{E}^{2}}=\frac{{b}_{0}D-{a}_{0}E}{{D}^{2}+{E}^{2}}\frac{{b}_{0}}{({a}_{0})^{2}+({b}_{0})^{2}},\nonumber\\ \end{eqnarray} \tag{ 44 }$$

or

$$\begin{eqnarray} &&\frac{{b}_{4}^{2}D-{a}_{4}{b}_{4}E}{({a}_{4})^{2}+({b}_{4})^{2}}=\frac{{b}_{0}^{2}D-{a}_{0}{b}_{0}E}{({a}_{0})^{2}+({b}_{0})^{2}}. \end{eqnarray} \tag{ 45 }$$

Substituting D = a₀a₄ − b₀b₄ + Ua₁ and E = a₀b₄ + a₄b₀ − Ub₁ one can easily show that the equality does not hold. Thus, our assumption is wrong and ${G}_{\sigma \sigma }^{\lt }(\omega )$ is not an imaginary function as it should be by definition.

If one is only interested in the Coulomb blockade regime, it is not necessary to go to the level of approximation presented here (which is essential to obtain the Kondo effect). For the Coulomb blockade regime one can turn to the approximation presented in [33], where on top of the approximations described above we also neglect the simultaneous hopping of electron pairs to and from the system. This approximation does not violet the symmetry relations of the single particle GF (see appendix A for further discussion), but as pointed above, it does not reproduce the Kondo peaks at low temperatures.

For the double Anderson model we follow the derivation given in [52], and define the following contour ordered GF

$$\begin{eqnarray} {G}_{\alpha \beta }^{\sigma \sigma }(\tau ,{\tau }^{\prime})&=&-\frac{\mathrm{i}}{\hbar }\langle {T}_{\mathrm{C}}{d}_{\alpha \sigma }(\tau ){d}_{\beta \sigma }^{\dagger }({\tau }^{\prime})\rangle , \end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray} {\mathbb{G}}_{\alpha \beta \gamma }^{s \sigma \sigma }(\tau ,{\tau }^{\prime})&=&-\frac{\mathrm{i}}{\hbar }\langle {T}_{\mathrm{C}}{n}_{\alpha s}(\tau ){d}_{\beta \sigma }(\tau ){d}_{\gamma \sigma }^{\dagger }({\tau }^{\prime})\rangle , \end{eqnarray} \tag{ 47 }$$

where $s=\sigma ,\bar {\sigma }$. The approximations used in [52] are: (a) neglect the simultaneous hopping of electron pairs to and from the system, (b) assume that ${F}_{2}(\tau ,{\tau }^{\prime})=-\frac{\mathrm{i}}{\hbar }\langle {T}_{\mathrm{C}}{c}_{k \sigma }(\tau )n(\tau ){d}_{\beta \sigma }^{\dagger }({\tau }^{\prime})\rangle \approx -\frac{\mathrm{i}}{\hbar }{\sum }_{\gamma =\alpha ,\beta }{t}_{k \gamma }^{\sigma }\int \mathrm{d}{\tau }_{1}{g}_{k}(\tau ,{\tau }_{1}) \langle {T}_{\mathrm{C}}{d}_{\gamma \sigma }({\tau }_{1})n({\tau }_{1}){d}_{\beta \sigma }^{\dagger }({\tau }^{\prime})\rangle$ where n(τ) is the number operator of one of the electrons of the system, and $(\mathrm{i}\hbar \frac{\partial }{\partial \tau }-{\varepsilon }_{k \sigma }){g}_{k}(\tau ,{\tau }_{1})=\delta (\tau -{\tau }_{1})$, and (c) higher order GFs of the form $-\frac{\mathrm{i}}{\hbar }\langle {T}_{\mathrm{C}}{n}_{\gamma \sigma }(\tau ){n}_{\delta \sigma }(\tau ){d}_{\alpha \sigma }(\tau ){d}_{\beta \sigma }^{\dagger }({\tau }^{\prime})\rangle$ are decoupled to $-\frac{\mathrm{i}}{\hbar }\langle {n}_{\gamma \sigma }(\tau )\rangle \langle {T}_{\mathrm{C}}{n}_{\delta \sigma }(\tau ){d}_{\alpha \sigma }(\tau ){d}_{\beta \sigma }^{\dagger }({\tau }^{\prime})\rangle \!-\!\frac{\mathrm{i}}{\hbar }\langle {n}_{\delta \sigma }(t)\rangle \langle {T}_{\mathrm{C}}{n}_{\gamma \sigma }(\tau ) {d}_{\alpha \sigma }(\tau ){d}_{\beta \sigma }^{\dagger }({\tau }^{\prime})\rangle$. These approximations lead to the following equations for the retarded and advanced GF (omitting (ω) for brevity):

$$\begin{eqnarray} &&\fl ({G}_{\alpha \beta }^{\sigma \sigma })^{r}=(\hbar \omega -{\varepsilon }_{\alpha \sigma }-{\Sigma }_{0}^{r})^{-1}\times [{\delta }_{\alpha \beta }+{h}_{\alpha \beta }^{\sigma }({G}_{\beta \beta }^{\sigma \sigma })^{r}\nonumber\\ &&+~{U}_{\alpha }({\mathbb{G}}_{\alpha \alpha \beta }^{\bar {\sigma }\sigma \sigma })^{r}+{V}_{\alpha \beta }^{\sigma \bar {\sigma }}({\mathbb{G}}_{\beta \alpha \beta }^{\bar {\sigma }\sigma \sigma })^{r}+{V}_{\alpha \beta }^{\sigma \sigma }({\mathbb{G}}_{\beta \alpha \beta }^{\sigma \sigma \sigma })^{r}], \end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray} &&\fl ({G}_{\beta \alpha }^{\sigma \sigma })^{a}=(\hbar \omega -{\varepsilon }_{\beta \sigma }-{\Sigma }_{0}^{a})^{-1}\times [{\delta }_{\beta \alpha }+{h}_{\beta \alpha }^{\sigma }({G}_{\alpha \alpha }^{\sigma \sigma })^{a}\nonumber\\ &&+~{U}_{\beta }({\mathbb{G}}_{\beta \beta \alpha }^{\bar {\sigma }\sigma \sigma })^{a}+{V}_{\beta \alpha }^{\sigma \bar {\sigma }}({\mathbb{G}}_{\alpha \beta \alpha }^{\bar {\sigma }\sigma \sigma })^{a}+{V}_{\beta \alpha }^{\sigma \sigma }({\mathbb{G}}_{\alpha \beta \alpha }^{\sigma \sigma \sigma })^{a}]. \end{eqnarray} \tag{ 49 }$$

The equations for the various $({\mathbb{G}}_{\alpha \beta \alpha }^{s \sigma \sigma }(\omega ))^{r,a}$ are provided in appendix B. Given the expressions for the retarded and advanced GF (equations (48) and (49)), one can show (see appendix B for a full proof) that $({G}_{\alpha \beta }^{\sigma \sigma }(\omega ))^{r}\not = (({G}_{\beta \alpha }^{\sigma \sigma }(\omega ))^{a})^{\ast }$. Moreover, we find that none of the symmetry relations in equation (18) hold. In the following section we propose a symmetrization scheme that restores all the symmetries of the single particle GF.

The customary route to calculate the NEGF is as follows: (a) calculate the retarded GF and use it to obtain the advanced GF (by demanding ${G}_{\alpha \beta }^{a}(\omega )=({G}_{\beta \alpha }^{r}(\omega ))^{\ast }$). (b) Calculate the lesser/greater GF and symmetrize the lesser/greater to fulfil the quantum Onsager relations [53], hence obeying ${G}_{\alpha \beta }^{\lt ,\gt }(\omega )=-({G}_{\beta \alpha }^{\lt ,\gt }(\omega ))^{\ast }$. In most applications of NEGF the advanced GF is not directly calculated and thus, the symmetry breakage does not always stand out. In fact, this common procedure restores the relation between the advanced and retarded GF and between the lesser/greater and their complex conjugate, but does not necessarily restore the relation ${G}_{\alpha \beta }^{r}(\omega )-{G}_{\alpha \beta }^{a}(\omega )={G}_{\alpha \beta }^{\gt }(\omega )-{G}_{\alpha \beta }^{\lt }(\omega )$. It can be shown that violation of the latter leads to violation of the fluctuation dissipation relation, G^< =− f_eq(ε − μ_eq)(G^r − G^a), at equilibrium. This oversimplified procedure can result in different values for the currents depending on how it is calculated, cf equation (1) or (2). It may also lead to finite currents at zero-bias voltage (see section 4.3 for more), which is physically incorrect.

In order to restore the symmetry relations that are imposed by the definitions of the GF (equation (18)), we suggest the following procedure:

• (i)  
  Calculate the retarded/advanced GFs matrices (G^r/G^a) separately and use them to define 'new' retarded/advanced GFs matrices ${\tilde {\mathbf{G}}}^{r}=0.5({\mathbf{G}}^{r}+({\mathbf{G}}^{a})^{\dagger })$ and ${\tilde {\mathbf{G}}}^{a}=0.5({\mathbf{G}}^{a}+({\mathbf{G}}^{r})^{\dagger })=({\tilde {\mathbf{G}}}^{r})^{\dagger }$.
• (ii)  
  Use the 'new' retarded/advanced GFs matrices $({\tilde {\mathbf{G}}}^{r}/{\tilde {\mathbf{G}}}^{a})$ to calculate the lesser/greater GFs matrices (G^</G^>). Again, use them to define 'new' lesser/greater GFs matrices ${\tilde {\mathbf{G}}}^{\lt ,\gt }=0.5({\mathbf{G}}^{\lt ,\gt }-({\mathbf{G}}^{\lt ,\gt })^{\dagger })$.
• (iii)  
  Calculate the two anti-Hermitian matrices $\mathbf{A}={\tilde {\mathbf{G}}}^{\gt }-{\tilde {\mathbf{G}}}^{\lt }$ and $\mathbf{B}={\tilde {\mathbf{G}}}^{r}-{\tilde {\mathbf{G}}}^{a}$. Define the difference anti-Hermitian matrix C = A − B, and redefine the retarded and advanced GFs ${\bar {\mathbf{G}}}^{r}={\tilde {\mathbf{G}}}^{r}+\mathbf{C}/2$, and ${\bar {\mathbf{G}}}^{a}={\tilde {\mathbf{G}}}^{a}-\mathbf{C}/2$.

The resulting GFs (${\bar {\mathbf{G}}}^{r},\hspace{0.167em} {\bar {\mathbf{G}}}^{a},\hspace{0.167em} {\tilde {\mathbf{G}}}^{\lt }$ and ${\tilde {\mathbf{G}}}^{\gt }$) obey all symmetry relations of equation (18) by construction. Note that if the original GFs obeyed the symmetry relations to begin with, our symmetrization procedure will not alter them in any way.

We now turn to perform detailed calculations for both the Anderson and double Anderson models. For the Anderson model, we use the closure described in section 3.2 while for the double Anderson model we use the closure described in section 3.3. The resulting EOMs were solved self-consistently in Fourier space with a frequency discretization of N_ω = 2¹⁴–2¹⁶ depending on the model parameters. Typically, <15 self-consistent iterations were needed to converge the results. Convergence was declared when the population values at subsequent iteration steps did not change within a predefined tolerance value chosen as 10⁻⁶. For each set of calculations we have applied the above symmetrization scheme and compared the results to those obtained without restoring symmetry, as detailed for each model.

First, we address the effects of symmetry breakage in the Anderson model. The closure used is sufficient to describe the appearance of the Kondo resonances at low temperatures, as seen in the upper panel of figure 2, where we plot the density of states as a function of energy for several temperatures, all calculated with symmetry restoration. The development of Kondo peaks in the density of states as the temperature decreases is clearly evident, signifying a regime of strong correlations which is qualitatively captured by the simple EOM approach when symmetry is restored.

Zoom In Zoom Out Reset image size

In the lower panel of figure 2 we show one of the main flaws of the EOM approach for the Anderson impurity model, where we plot the value of 〈n_↑〉 as a function of the source–drain bias voltage with and without symmetry restoration. The most notable effect is the appearance of an imaginary portion to 〈n_↑〉 as the source–drain bias voltage is increased. To obtain the results, without symmetry restoration, only the real part of 〈n_σ〉 was used to converge the self-consistent equations for the GFs. By applying the symmetrization scheme proposed in section 4.1 to the lesser GF calculated in section 3.2, we restore the relation ${G}_{\alpha \beta }^{\lt ,\gt }(\omega )=-({G}_{\beta \alpha }^{\lt ,\gt }(\omega ))^{\ast }$. This is sufficient to obtain a real value for 〈n_↑〉, as clearly shown in the lower panel of figure 2. All other symmetry relation are not violated here and thus, our symmetrization procedure does not affect them at all. Interestingly, taking only the real part of 〈n_σ〉 provides identical results when compared to the results obtained after the full symmetrization procedure. However, this is only true for the simple case of the single site impurity model and does not hold for more complex systems.

We now turn to discuss the impact of symmetry breaking for the double Anderson model. This system is more involved compared to the single site Anderson model and thus, the level of closure used is somewhat simpler, as explained in section 3.3. While for the case of a single site Anderson model only the relation ${G}_{\alpha \beta }^{\lt ,\gt }(\omega )=-({G}_{\beta \alpha }^{\lt ,\gt }(\omega ))^{\ast }$ breaks down, in the double Anderson model we find that all 3 symmetries described by equation (18) are violated. This can be traced to the more complex form of the Hamiltonian for the double Anderson model, where each site is only coupled to one of the leads and transport in enabled by the direct hopping term between the two sites.

Similar to the case of the Anderson model, as a result of symmetry breaking the occupation of the levels 〈n_ασ〉 is a complex number. In addition, the coherences, ${\rho }_{\alpha \beta }^{\sigma \sigma }=-\frac{\mathrm{i}\hbar }{2 \pi }\int \nolimits _{-\infty }^{\infty }({G}_{\alpha \beta }^{\sigma \sigma }(\omega ))^{\lt }\hspace{0.167em} \mathrm{d}\omega$, should also fulfil certain symmetry relations, such as ${\rho }_{\alpha \beta }^{\sigma \sigma }=({\rho }_{\beta \alpha }^{\sigma \sigma })^{\ast }$. In figure 3 we plot the real and imaginary parts of ${\rho }_{\alpha \beta }^{\sigma \sigma }$ and ${\rho }_{\beta \alpha }^{\sigma \sigma }$ for the case where the symmetry procedure has been applied (left panels) and for the bare case (right panels). The upper panels show the imaginary part of ${\rho }_{\alpha \beta }^{\sigma \sigma }$ and ${\rho }_{\beta \alpha }^{\sigma \sigma }$, which should show a mirror reflection about the zero axis (shown as thin solid line). This is, indeed, the case when symmetry is restored, however, it is destroyed when symmetry breaks down, in particular as the source–drain bias increases. A more dramatic effect is shown for the real part of ${\rho }_{\alpha \beta }^{\sigma \sigma }$ and ${\rho }_{\beta \alpha }^{\sigma \sigma }$ (lower panels). The two curves representing $\mathrm{Re}({\rho }_{\alpha \beta }^{\sigma \sigma })$ and $\mathrm{Re}({\rho }_{\beta \alpha }^{\sigma \sigma })$ should be identical (left panel when symmetry is restored) but are quite distinct when symmetry is not obeyed (right panel).

Zoom In Zoom Out Reset image size

In figure 4 we plot the current as a function of the source–drain bias voltage for the double Anderson model. The current can be obtained from equation (1) (dashed line) or from equation (2) (dotted curve). In the limit of infinite hierarchy in the EOM approach the two formulae should coincide. However, when approximations are introduced or when the hierarchy is truncated, the calculation of the current based on the two different formulae will coincide only if the symmetry relation $({G}_{\alpha \beta }^{\sigma \tau }(\omega ))^{r}-({G}_{\alpha \beta }^{\sigma \tau }(\omega ))^{a}=({G}_{\alpha \beta }^{\sigma \tau }(\omega ))^{\gt }-({G}_{\alpha \beta }^{\sigma \tau }(\omega ))^{\lt }$ is preserved. Indeed, in the case of a single site Anderson model, even if symmetry is not restored, this relation holds and the two calculations yield identical values for the current. However, in the present case, all 3 symmetry relations are broken and thus, equations (1) and (2) give different results for the current, as clearly evident in figure 4. More significantly is the fact that equation (1) produces a finite value for the current even when the bias is zero, indicating the break down of the fluctuation dissipation relation. When symmetry is restored (solid curve) the two calculations are identical, as they should be, and the violation of the fluctuation dissipation relation is also resolved.

Zoom In Zoom Out Reset image size

The symmetrization scheme proposed here is not a 'magic cure' and, in fact, does not resolve all issues of mater. It is well known that the lesser and greater GFs should obey a simple sum rule where the integral over the difference of their diagonal elements should always sum to 1:

$$\begin{eqnarray} &&{S}_{\alpha \sigma }=\frac{-\mathrm{i}\hbar }{2 \pi }\int \mathrm{d}\varepsilon (({G}_{\alpha \alpha }^{\sigma \sigma }(\epsilon ))^{\lt }-({G}_{\alpha \alpha }^{\sigma \sigma }(\epsilon ))^{\gt })=1. \end{eqnarray} \tag{ 50 }$$

In figure 5 we plot the sum rule as given by equation (50) for the double Anderson model where symmetry has been restored. A similar plot for the single site Anderson model yields a value of 1 regardless of whether symmetry has been restored or not within the closure discussed above. However, in the case of the more evolved double Anderson model, even when symmetry is restored and the GFs obey all 3 relations described in equation (18), the sum rule is violated. Nonetheless, the sum ∑_ασS_ασ = N_e, where N_e is the total number of electrons in the system at maximal occupancy, is indeed preserved when symmetrization is restored.

Zoom In Zoom Out Reset image size