Many-body current formula and current conservation for non-equilibrium fully interacting nanojunctions

The Green function on the Keldysh contour C_K is defined as

$$\begin{equation} G(1,2) = - {\rm i} \langle \mathcal {T}_{C_K} \Psi (1) \Psi ^\dag (2) \rangle , \end{equation} \tag{ 1 }$$

where (1, 2) stands for a composite index for space $\mathbf {x}_{1,2}$ (i.e. states λ, n or ρ in the L, C or R regions, respectively) and time τ_{1, 2} on the time-loop contour C_K. The time ordering $\mathcal {T}_{C_K}$ of the product of fermion creation (Ψ†) and annihilation (Ψ) quantum fields is performed on the time-loop contour C_K [36–39].

The Green function obeys the equation of motion on the contour C_K [37, 39]:

$$\begin{equation} [{\rm i}\partial _{\tau _1} - \bar{h}_0(1)]G(1,1^\prime ) = \delta (1 - 1^\prime ) + \int {\rm d}3 \Sigma (1,3)G(3,1^\prime ), \end{equation} \tag{ 2 }$$

and the adjoint equation reads

$$\begin{equation} [-{\rm i}\partial _{\tau _1^\prime } - \bar{h}_0(1^\prime )]G(1,1^\prime ) = \delta (1 - 1^\prime ) + \int {\rm d}2 G(1,2) \Sigma (2,1^\prime ), \end{equation} \tag{ 3 }$$

where $\bar{h}_0(1)$ is the non-interacting Hamiltonian.

Using the continuity equation $\nabla \vec{j}(1) + \partial _t n(1) = 0$, one can write down the current through the interface between the L and C regions. It should be noted that considering an interface between the L (R) and C regions is merely a virtual partitioning for mathematical convenience—the interactions are present throughout the entire system. It is also useful to consider such LC/RC interfaces in order to connect our results to previously derived expressions within the partitioning scheme. The location of these interfaces is also arbitrary, but for convenience in future numerical computation, different choices are possible: for example, at the contacts between the leads and the ends of the molecule, or at the contacts between the leads and the so-called extended molecule, which already contains part of the leads.

After integration of the continuity equation over the (half) space of the L region, one obtains the following expression for the current flowing through the left interface (from L to C):

$$\begin{equation} I_L(t) = - e \frac{\rm d}{{\rm d}t} \langle \hat{N}_L(t) \rangle , \end{equation} \tag{ 4 }$$

where $\langle \hat{N}_L(t) \rangle$ is the total number of electrons in the L region. It is obtained from the lesser Green function as

$$\begin{equation} \langle \hat{N}_L(t) \rangle = \sum _\lambda -{\rm i} G^<_{\lambda \lambda }(t,t). \end{equation} \tag{ 5 }$$

To calculate the time derivative of G^<_λλ(t, t), we first go back to the full two-time dependence of G^<, and then take the equal-time limit after performing the derivatives:

$$\begin{equation} \frac{\rm d}{{\rm d}t}\langle \hat{N}_L(t) \rangle = \sum _\lambda \left( -{\rm i} \frac{\rm d}{{\rm d}t_1} G^<_{\lambda \lambda }(t_1,t) - {\rm i} \frac{d}{{\rm d}t_2} G^<_{\lambda \lambda }(t,t_2) \right)_{t_1=t_2= t}. \end{equation} \tag{ 6 }$$

Using the equations of motion, equations (2) and (3), for the Green functions on C_K, one obtains the current I_L as

$$\begin{equation} I_L(t) = \frac{e}{\hbar } {\rm Tr}_\lambda [ (\Sigma G)^<(t,t) - (G \Sigma )^<(t,t) ], \end{equation} \tag{ 7 }$$

where we have re-introduced the ℏ to have the correct units of current and conductance.

Using the rules of analytical continuation (see appendix B), we get

$$\begin{eqnarray} &&\fl {\rm Tr}_\lambda [\ldots ] = {\rm Tr}_\lambda [ (\Sigma ^{{\rm MB},<} G^a) (t,t) + ( ( V_{LC} + \Sigma ^{{\rm MB},r})(t,t) G^< ) \nonumber\\ &&-\, ( G^< (V_{LC} + \Sigma ^{{\rm MB},a} ) )(t,t) - (G^r\Sigma ^{{\rm MB},<})(t,t) ] . \end{eqnarray} \tag{ 8 }$$

There are no lesser and greater components for V_LC since its time dependence is local, i.e. V_LC(t, t') = V_LC(t)δ(t − t').

In the steady state, all double-time quantities X(t, t') depend only on the time difference X(t − t'). One obtains the following expression for the current I_L after the Fourier transform

$$\begin{equation} I_L = \frac{e}{\hbar } \int \frac{{\rm d}\omega }{2\pi } {\rm Tr}_\lambda [ (\Sigma (\omega ) G(\omega ))^< - (G(\omega ) \Sigma (\omega ))^< ] \end{equation} \tag{ 9 }$$

or equivalently

$$\begin{eqnarray} &&\fl I_L = \frac{e}{\hbar } \int \frac{{\rm d}\omega }{2\pi } {\rm Tr}_\lambda [ V_{LC}G^<(\omega ) - G^<(\omega ) V_{LC} + \Sigma ^{{\rm MB},<}(\omega ) G^a(\omega ) + \Sigma ^{{\rm MB},r}(\omega ) G^<(\omega ) \nonumber\\ &&-\, G^<(\omega )\Sigma ^{{\rm MB},a}(\omega ) - G^r(\omega ) \Sigma ^{{\rm MB},<}(\omega ) ]. \end{eqnarray} \tag{ 10 }$$

Now that we have sorted out the Keldysh components, we have to sort out the index (matrix elements) parts. First we concentrate on Tr_λ[(Σ^MBG)^<] and Tr_λ[(GΣ^MB)^<]. We thus have (we use the symbol $\bullet \hspace{-8.25128pt}\sum$ for summation to have a better graphical distinction between the sums and the self-energies Σ):

$$\begin{eqnarray} &&{\rm Tr}_\lambda [ (\Sigma ^{{\rm MB}} G)^< ] = \bullet \hspace{-9.95845pt}\sum _{\lambda ,n} \big(\Sigma ^{{\rm MB}}_{\lambda n} G_{n \lambda }\big)^< + \bullet \hspace{-9.95845pt}\sum _{\lambda ,\lambda ^{\prime }} \big(\Sigma ^{{\rm MB}}_{\lambda \lambda ^{\prime }} G_{\lambda ^{\prime } \lambda }\big)^< + \bullet \hspace{-9.95845pt}\sum _{\lambda ,\rho } \big(\Sigma ^{{\rm MB}}_{\lambda \rho } G_{\rho \lambda }\big)^<, \end{eqnarray} \tag{ 11 }$$

and similarly for Tr_λ[(GΣ^MB)^<]:

$$\begin{eqnarray} &&{\rm Tr}_\lambda [ (G \Sigma ^{{\rm MB}})^< ] = \bullet \hspace{-9.95845pt}\sum _{\lambda ,n} \big(G_{\lambda n} \Sigma ^{{\rm MB}}_{n \lambda }\big)^< + \bullet \hspace{-9.95845pt}\sum _{\lambda ,\lambda ^{\prime }} \big(G_{\lambda \lambda ^{\prime }} \Sigma ^{{\rm MB}}_{\lambda ^{\prime } \lambda }\big)^< + \bullet \hspace{-9.95845pt}\sum _{\lambda ,\rho } \big(G_{\lambda \rho } \Sigma ^{{\rm MB}}_{\rho \lambda }\big)^<. \end{eqnarray} \tag{ 12 }$$

In this present version of the theory, we assume that the matrix elements Σ^MB_ρλ and Σ^MB_λρ vanish, as there is no direct interaction between the left and right electrode. Because of the geometry and the heterogeneity of the nanojunctions and because of the different dimensionality of the leads and the central region, we assume that there is an effective screening of the Coulomb interaction so that the electrons of the L and R leads do not interact directly via any Σ^MB_ρλ or Σ^MB_λρ matrix elements. The distance $\vert \mathbf {x}_\lambda -\mathbf {x}_\rho \vert$ between two points in the L and R leads respectively is large enough so that the spatial decay of $\Sigma ^{\rm MB}_{\lambda \rho } = \Sigma ^{\rm MB}(\vert \mathbf {x}_\lambda -\mathbf {x}_\rho \vert )$ make the contribution of Σ^MB_λρ zero or negligible. In other words, the presence of the L and R electrodes affect directly the central region C via Σ^MB_LC and Σ^MB_RC. However, the L electrode do not affect directly the R electrode (and vice versa), but only indirectly via exchange and correlation effects involving the states of the central region C. This assumption seems to be valid when the size of the central region is of the order of several atoms (i.e. a molecule). If the central region were to be a single atomic impurity coupled to two continuum of the delocalized electron state, this assumption would not be valid.

We now turn to the evaluation of the sums in equations (11) and (12). First we concentrate on the $\bullet \hspace{-8.25128pt}\sum _{\lambda ,\lambda ^{\prime }}$ sums

$$\begin{eqnarray} &&\fl \sum _{\lambda ,\lambda ^{\prime }} \big(\Sigma ^{{\rm MB}}_{\lambda \lambda ^{\prime }} G_{\lambda ^{\prime } \lambda }\big)^< - \bullet \hspace{-9.95845pt}\sum _{\lambda ,\lambda ^{\prime }} \big(G_{\lambda \lambda ^{\prime }} \Sigma ^{{\rm MB}}_{\lambda ^{\prime } \lambda }\big)^< \nonumber\\ &&\ms = \bullet \hspace{-9.95845pt}\sum _{\lambda ,\lambda ^{\prime }} \Sigma ^{{\rm MB}, <}_{\lambda \lambda ^{\prime }} G_{\lambda ^{\prime } \lambda }^a + \Sigma ^{{\rm MB},r}_{\lambda \lambda ^{\prime }} G_{\lambda ^{\prime } \lambda }^< - G^<_{\lambda \lambda ^{\prime }} \Sigma ^{{\rm MB},a}_{\lambda ^{\prime } \lambda } - G^r_{\lambda \lambda ^{\prime }} \Sigma ^{{\rm MB},<}_{\lambda ^{\prime } \lambda } \nonumber\\ &&\ms = \bullet \hspace{-9.95845pt}\sum _{\lambda ,\lambda ^{\prime }} \Sigma ^{{\rm MB}, <}_{\lambda \lambda ^{\prime }} \big( G_{\lambda ^{\prime } \lambda }^a - G_{\lambda ^{\prime } \lambda }^r \big) + \big( \Sigma ^{{\rm MB},r}_{\lambda \lambda ^{\prime }} - \Sigma ^{{\rm MB},a}_{\lambda ^{\prime } \lambda }\big) G_{\lambda ^{\prime } \lambda }^< \nonumber\\ &&\ms = \bullet \hspace{-9.95845pt}\sum _{\lambda ,\lambda ^{\prime }} \Sigma ^{{\rm MB}, <}_{\lambda \lambda ^{\prime }} (G^<-G^>)_{\lambda ^{\prime } \lambda } + ( \Sigma ^{{\rm MB},>} - \Sigma ^{{\rm MB},<} )_{\lambda \lambda ^{\prime }} G_{\lambda ^{\prime } \lambda }^< \nonumber\\ &&\ms = \bullet \hspace{-9.95845pt}\sum _{\lambda ,\lambda ^{\prime }} \Sigma ^{{\rm MB}, >}_{\lambda \lambda ^{\prime }} G^<_{\lambda ^{\prime } \lambda } - \Sigma ^{{\rm MB},<}_{\lambda \lambda ^{\prime }} G_{\lambda ^{\prime } \lambda }^> = {\rm Tr}_{\lambda } \left[ \Sigma ^{{\rm MB}, >}_L G^<_L - \Sigma ^{{\rm MB},<}_L G^>_L \right]. \end{eqnarray} \tag{ 13 }$$

In the second line, we have used the rules of analytical continuation. In the third, we have used the equivalent of cyclic permutation in the calculation of a trace, i.e. swapping the indices λ and λ' in the last two terms. This is possible here since the sums and all matrix elements are defined in the single subspace of the L electrode. The final result looks like the collision terms usually obtained in the derivation of the generalized Boltzmann equation from quantum kinetic theory. They correspond to the particle production (scattering-in) and absorption or hole production (scattering-out) related to inelastic processes (i.e. non-diagonal elements of the self-energy on the time-loop contour Σ^<) occurring in the left electrode. It has also been shown that the integration of such term vanishes as a result of the gauge degree of freedom [40]. In section 3.4, we use another route to show how and why these terms can vanish by using generalized non-equilibrium distribution functions.

Now let us concentrate on the $\bullet \hspace{-8.25128pt}\sum _{\lambda ,n}$ sums in equation (11) and equation (12). Once more using the rules of analytical continuation, we can see that we need the knowledge of the following Green's functions matrix elements: G^<_nλ, G^<_λn, G^a_nλ and G^r_λn. For this we use the Dyson-like equation defined for the non-diagonal elements: G^<_nλ = 〈n|(GΣg)^<|λ〉. We will not go into the detail of the full calculations for all four Green's function matrix elements; instead we concentrate on one matrix element 〈n|(GΣg)^<|λ〉 to show the mechanism of the derivation

$$\begin{eqnarray} &&\fl G^<_{n\lambda } = \langle n\vert (G\Sigma g)^<\vert \lambda \rangle =\langle n\vert G^r \Sigma ^< g^r + G^< \Sigma ^a g^a + G^r \Sigma ^r g^< \vert \lambda \rangle \nonumber\\ &&=\bullet \hspace{-18.25128pt}\sum _{\lambda _1,\lambda _2,\rho ,m} G^r_{n\lambda _1} \Sigma ^<_{\lambda _1 \lambda _2} g^a_{\lambda _2 \lambda } + G^r_{n\rho } \Sigma ^<_{\rho \lambda _2} g^a_{\lambda _2 \lambda } + G^r_{nm} \Sigma ^<_{m \lambda _2} g^a_{\lambda _2 \lambda } \nonumber\\ &&+\, G^<_{n\lambda _1} \Sigma ^a_{\lambda _1 \lambda _2} g^a_{\lambda _2 \lambda } + G^<_{n\rho } \Sigma ^a_{\rho \lambda _2} g^a_{\lambda _2 \lambda } + G^<_{nm} \Sigma ^a_{m \lambda _2} g^a_{\lambda _2 \lambda } \nonumber\\ &&+\, G^r_{n\lambda _1}\Sigma ^r_{\lambda _1 \lambda _2} g^<_{\lambda _2 \lambda } + G^r_{n\rho } \Sigma ^r_{\rho \lambda _2} g^<_{\lambda _2 \lambda } + G^r_{nm}\ \Sigma ^r_{m \lambda _2}\ g^<_{\lambda _2 \lambda } , \end{eqnarray} \tag{ 14 }$$

with Σ^a/r_mλ = V_mλ + Σ^{MB, a/r}_mλ and Σ^<_mλ = Σ^{MB, <}_mλ. One has to keep in mind that we use a model such that Σ^x_ρλ = 0.

As we show in detail in section 4.1, the interaction defined within the subspace of the L lead $\Sigma ^{a/r}_{\lambda _1 \lambda _2}$ can be factorized out and included in the renormalized Green functions $\tilde{g}^{a/r,<}_{\lambda _1 \lambda _2}$ of the L lead. Hence, the matrix elements G^<_nλ = 〈n|(GΣg)^<|λ〉 can be recast as $G^<_{n\lambda }=\langle n\vert (G_C\ \Sigma _{CL}\ \tilde{g}_L)^<\vert \lambda \rangle$ with

$$\begin{equation} G^<_{n\lambda }= \bullet \hspace{-9.95845pt}\sum _{m,\lambda ^{\prime }} G^r_{nm}\ \Sigma ^r_{m \lambda ^{\prime }}\ \tilde{g}^<_{\lambda ^{\prime } \lambda } + G^r_{nm}\ \Sigma ^{{\rm MB},<}_{m \lambda ^{\prime }}\ \tilde{g}^a_{\lambda ^{\prime } \lambda } + G^<_{nm}\ \Sigma ^a_{m \lambda _2}\ \tilde{g}^a_{\lambda ^{\prime } \lambda } . \end{equation} \tag{ 15 }$$

Similarly, we find that

$$\begin{eqnarray} &&G^<_{\lambda n} = \langle \lambda \vert (\tilde{g}_L\ \Sigma _{LC}\ G_C)^<\vert n \rangle , \nonumber\\ &&G^a_{n\lambda } = \langle n\vert (G_C\ \Sigma _{CL}\ \tilde{g}_L)^a\vert \lambda \rangle , \\ &&G^r_{\lambda n} = \langle \lambda \vert (\tilde{g}_L\ \Sigma _{LC}\ G_C)^r\vert n \rangle . \nonumber \end{eqnarray} \tag{ 16 }$$

Finally by using the rules of analytical continuation for the above matrix elements, we find that the current flowing through the left L interface (10) can be recast as

$$\begin{eqnarray} &&\fl I_L \,{=}\, \frac{e}{\hbar } \!\int \frac{{\rm d}\omega }{2\pi } {\rm Tr}_n\big[ G^r_C \tilde{\Upsilon }^{L,l}_C + G^a_C \big(\tilde{\Upsilon }^{L,l}_C\big)^\dagger + G^<_C \big( \tilde{\Upsilon }^L_C - \big(\tilde{\Upsilon }^L_C\big)^\dagger \big) \big] + {\rm Tr}_\lambda \left[\Sigma ^{{\rm MB},>}_{L} G^<_{L} - \Sigma ^{{\rm MB},<}_{L} G^>_{L} \right], \nonumber\\ \end{eqnarray} \tag{ 17 }$$

where

$$\begin{eqnarray} \tilde{\Upsilon }^L_C &=& \Sigma ^a_{CL}\ \tilde{g}^a_L\ \Sigma ^r_{LC} , \nonumber\\ \tilde{\Upsilon }^{L,l}_C &=& \Sigma ^<_{CL} \left( \tilde{g}^a_L - \tilde{g}^r_L \right) \Sigma ^r_{LC} + \Sigma ^r_{CL}\ \tilde{g}^<_L\ \Sigma ^r_{LC} \nonumber\\ &=& (\Sigma \tilde{g})^<_{CL}\ \Sigma ^r_{LC} - \Sigma ^<_{CL}\ (\tilde{g} \Sigma )^r_{LC} \end{eqnarray} \tag{ 18 }$$

and

$$\begin{eqnarray} \big(\tilde{\Upsilon }^L_R\big)^\dag &=& \Sigma ^a_{CL}\ \tilde{g}^r_L\ \Sigma ^r_{LC} , \nonumber\\ \big(\tilde{\Upsilon }^{L,l}_C\big)^\dag &=& \Sigma ^a_{CL} \left( \tilde{g}^a_L - \tilde{g}^r_L \right) \Sigma ^<_{LC} - \Sigma ^a_{CL}\ \tilde{g}^<_L\ \Sigma ^a_{LC} \nonumber\\ &=& (\Sigma \tilde{g})^a_{CL}\ \Sigma ^<_{LC} - \Sigma ^a_{CL}\ (\tilde{g} \Sigma )^<_{LC} . \end{eqnarray} \tag{ 19 }$$

Equations (17), (18) and (19) are the main results of this work for the current formula. The current I_L flowing at the left LC interface is given by two traces: the first trace is a generalization of the Meir and Wingreen expression of the current to the cases where there are both interactions within the left electrode and crossing at the LC contact. The second trace is a term related to inelastic transport effects involving summation over the left electrode states/sites only.

An expression similar to equation (17) can be obtained for the current I_R flowing at the right RC interface, by swapping the index L↔R and changing the sign to keep the same convention for positive current flowing from the left to right direction. We find

$$\begin{eqnarray} &&\fl I_R = - \frac{e}{\hbar } \int \frac{{\rm d}\omega }{2\pi } {\rm Tr}_n\big[ G^r_C \tilde{\Upsilon }^{R,l}_C + G^a_C \big(\tilde{\Upsilon }^{R,l}_C\big)^\dagger + G^<_C \big( \tilde{\Upsilon }^R_C - (\tilde{\Upsilon }^R_C)^\dagger \big) \big] \nonumber\\ &&+\, {\rm Tr}_\rho \left[\Sigma ^{{\rm MB},>}_{R} G^<_{R} - \Sigma ^{{\rm MB},<}_{R} G^>_{R} \right]. \end{eqnarray} \tag{ 20 }$$

We now comment on the physical meaning and implication of the new ϒ^α_C quantities and the different traces entering the definition of the current.

From the definition of $\tilde{\Upsilon }^{L,l}_C$, equation (18), we can also define the quantity $\tilde{\Upsilon }^{L,g}_C$ using the greater components Σ^>_CL and $\tilde{g}^r_L$ such as

$$\begin{equation} \tilde{\Upsilon }^{L,g}_{C} = \Sigma ^>_{CL} \left( \tilde{g}^a_L - \tilde{g}^r_L \right) \Sigma ^r_{LC} + \Sigma ^r_{CL}\ \tilde{g}^>_L\ \Sigma ^r_{LC} . \end{equation} \tag{ 21 }$$

It is now easy to show that $\tilde{\Upsilon }^l_{LC}$ and $\tilde{\Upsilon }^g_{LC}$ are related to each other by

$$\begin{equation} \tilde{\Upsilon }^{L,g}_{C} - \tilde{\Upsilon }^{L,l}_{C} = \big(\tilde{\Upsilon }^L_C\big)^\dag - \tilde{\Upsilon }^L_C . \end{equation} \tag{ 22 }$$

This is a very interesting relationship in the sense that the $\tilde{\Upsilon }^L_C$ quantities obey a relation of the type (greater) − (lesser) = (retarded) − (advanced) as for conventional self-energies or Green's functions, though $\tilde{\Upsilon }^{L,l}_C$ and $\tilde{\Upsilon }^{L,g}_C$ are not proper lesser and greater quantities. It should also be noted that the different $\tilde{\Upsilon }^L_C(\omega )$ play a similar role as the lead self-energies (or embedding potentials [41–43]), defined as Σ^{L, x}_C(ω) = V_CL g^x_L(ω) V_LC when the interactions are not crossing at the contacts. However, they are not simply related to the straightforward generalization of these embedding potentials. The latter have the following form (see also the expression for the Green functions given in section 4)

$$\begin{equation} \tilde{Y}^{L,x}_C(\omega ) = (\Sigma _{CL}(\omega )\ \tilde{g}_L(\omega )\ \Sigma _{LC}(\omega ))^x \end{equation} \tag{ 23 }$$

with Σ_LC/CL(ω) = V_LC/CL + Σ^MB_LC/CL(ω) and x = ( >, <, r, a). The rules of analytical continuation for $\tilde{Y}^{L,<}_C$ and $\tilde{Y}^{L,r}_C$ do not give the same expression as for $\tilde{\Upsilon }^{L,l}_C$ or $\tilde{\Upsilon }^L_C$. For example, $\tilde{Y}^{L,r}_C= \Sigma ^r_{CL}\ \tilde{g}^r_L\ \Sigma ^r_{LC} \ne \tilde{\Upsilon }^L_C$ and $\tilde{Y}^{L,<}_C \ne \tilde{\Upsilon }^{L,l}_C$.

This leads us to an important proof of this work: because of (a) the very existence of the interaction crossing at the contact, of (b) the fact that Σ^a_Lα/αC ≠ Σ^r_Lα/αC (as opposed to the non-interacting case where V^a_Lα/αC = V^r_Lα/αC = V_Lα/αC, and of (c) the rules of analytical continuation for products of three quantities, the usual cyclic permutation used in the calculation of the trace Tr_λ[(ΣG)^< − (GΣ)^<] cannot be used to transform the initial trace over {λ} onto a trace over {n} only. Therefore, the current I_L at the LC contact cannot be expressed simply in terms of the generalized embedding potentials $\tilde{Y}^{L,x}_C$; hence, the introduction of the $\tilde{\Upsilon }^L_C$ quantities. The current is not obtained from a straightforward generalization of the Meir and Wingreen formula using the embedding potentials $\tilde{Y}^{L,x}_C$:

$$\begin{equation} I_L \ne \frac{e}{\hbar } \int \frac{{\rm d}\omega }{2\pi } {\rm Tr}_n \big[Y^{L,<}_C G^>_C - Y^{L,>}_C G^<_C\big] + {\rm Tr}_\lambda \big[\Sigma ^{{\rm MB}>}_L G^<_L - \Sigma ^{{\rm MB}<}_L G^>_L\big]. \end{equation} \tag{ 24 }$$

Now we consider the terms Tr_λ[Σ^{MB, >}_LG^<_L − Σ^{MB, <}_LG^>_L] and show in which conditions this trace vanishes. For this we introduce the non-equilibrium distribution functions f^<_L(ω) and f^{MB, <}_L(ω) defined from the generalized Kadanoff–Baym ansatz [44] as follows:

$$\begin{eqnarray} &&G^<_L(\omega ) = f^<_L(\omega ) G^a_L(\omega ) - G^r_L(\omega ) f^<_L(\omega ) \\ &&\ms \Sigma ^{{\rm MB},<}_L(\omega ) = f^{{\rm int},<}_L(\omega ) \Sigma ^{{\rm MB},a}_L(\omega ) - \Sigma ^{{\rm MB},r}_L(\omega ) f^{{\rm int},<}_L(\omega ) , \end{eqnarray} \tag{ 25 }$$

and similarly for f^>_L and f^{MB, >}_L, obtained from G^>_L and Σ^{MB, >}_L. These distribution functions follow the conditions f^>_L = f^<_L − 1_L and f^{int, <}_L = f^{int, <}_L − 1_L (where the identity matrix in the L region is $1_L=\delta _{\lambda \lambda ^{\prime }}$ so that the usual relationships between the different Green functions (and self-energies ) X^r − X^a = X^> − X^< still hold.

For convenience, we consider for the moment that the non-equilibrium distribution functions are diagonal in the corresponding subspace, i.e. $f^<_{\lambda \lambda ^{\prime }} = f^<_\lambda \delta _{\lambda \lambda ^{\prime }}$. However, there is no formal difficulty to deal with a full, non-diagonal, dependence of these density-matrix-like distribution functions.

Introducing the definition of the non-equilibrium distribution functions in equation (13), we end up, after lengthy (but rather trivial) calculations, with

$$\begin{equation} \fl {\rm Tr}_{\lambda } \left[ \Sigma ^{{\rm MB}, >}_L G^<_L - \Sigma ^{{\rm MB},<}_L G^>_L \right] = (2\pi )^2 {\rm Tr}_\lambda \big[\big(f^<_L - f^{{\rm int},<}_L\big) A^\Sigma _L(\omega ) A^G_L(\omega ) \big], \end{equation} \tag{ 26 }$$

where the respective spectral functions are obtained from

$$\begin{equation} 2\pi {\rm i} A^X_L(\omega ) = X^a_L(\omega ) - X^r_L(\omega ), \end{equation} \tag{ 27 }$$

with X(ω) = G(ω) or Σ(ω).

Equation (26) shows that the collision term is a measure of the deviation between the two distribution functions f^<_L(ω) and f^{int, <}_L(ω). At equilibrium, because all distribution functions are equal to the Fermi distribution f^<_L = f^{int, <}_L = f^eq, this collision term vanishes. At non-equilibrium this is not generally the case.

One can now imagine the following case: the indices λ represent a spatial location or a localized electronic state on a lattice point in the left electrode. For λ located well inside the electrode, the system is in its local equilibrium (the Fermi distribution with a Fermi level shifted by the left bias) and both distribution functions f^<_L and f^{int, <}_L are equal to the left Fermi distribution; hence, these indices do not contribute to the trace. Only the lattice points (or states) that are not in local quasi-equilibrium, i.e. those close enough to the central region to experience the potential drop and the interaction effects at the contact and beyond will contribute to the trace. From a computational point of view, this is good news in the sense that one is not obliged to perform the summation in Tr_λ[...] over all the infinite λ indices of the semi-infinite left electrode. Only the states/sites for which f^<_λ − f^{int, <}_λ ≠ 0 will contribute to the trace. So if we choose the location of the LC interface far enough inside the left electrode of the real system, then Tr_λ[...] = 0. In this sense, we have just provided a formal proof of the concept of the so-called extended molecule [45, 10, 9, 6, 46, 47]. The extended molecule not only represents the central region C and consists of the molecule itself but also a part of the left and right electrodes to which the molecule is connected. The concept of the extended molecule has been introduced empirically in realistic calculations of molecular junctions in order to avoid any problems with the asymptotic behaviour of the electrostatic potential in the bias of finite (not small) applied bias.

The value of Tr_λ[...] can also be understood as a measure of the 'efficiency' of the location of the LC interface in the L electrode. The larger (smaller) the value is, the farther (closer) from local equilibrium the LC interface is. This measure can help in finding a good compromise between having a sufficiently large extended molecule that is nonetheless small enough for tractable numerical calculations.

Using different approximations (for example, single-particle approaches, mean-field theories, interactions localized in the central region only) for the different interacting self-energies, we can recover from equation (17) all the previously derived current expressions for non-equilibrium nanojunctions. We have analysed all these connections in detail in [35] and we will not repeat the analysis here.

However, as will be useful below, we now briefly recall that equation (17) bears some resemblance to the current expression derived by Meir and Wingreen [25] in the case of interaction present in the central region only:

$$\begin{equation} I_L^{\rm MW} = \frac{{\rm i} e}{\hbar } \int \frac{{\rm d}\omega }{2\pi }\ {\rm Tr}_{n} \left[ f_L \big(G^r_C - G^a_C\big) \Gamma _C^L + G^<_C \Gamma _C^L \right] . \end{equation} \tag{ 28 }$$

The connection becomes more apparent when we use the definitions if_LΓ^L_C = V_CL g^<_L V_LC = Σ^{L, <}_C = −(Σ^{L, <}_C)† and iΓ^L_C = V_CL(g^a_L − g^r_L)V_LC = Σ^{L, a}_C − Σ^{L, r}_C. Hence, I^MW_L becomes

$$\begin{eqnarray} &&I_L^{\rm MW} = \frac{e}{\hbar } \int \frac{{\rm d}\omega }{2\pi } {\rm Tr}_{n} \big[ G^r_C \Sigma ^{L,<}_C + G^a_C \big(\Sigma ^{L,<}_C\big)^\dag + G^<_C \big(\Sigma ^{L,a}_C - \Sigma ^{L,r}_C\big) \big] . \end{eqnarray} \tag{ 29 }$$

One can see by comparing equations (17) and (29) that the quantities $\tilde{\Upsilon }^L_C$, $\big(\tilde{\Upsilon }^L_C\big)^\dag$ and $\tilde{\Upsilon }^{L,l}_C$ play the role of the L lead self-energies Σ^{L, a}_C, Σ^{L, r}_C and Σ^{L, <}_C, respectively, in the cases where the interactions cross at the L interface. However, as we have discussed above, the quantities $\tilde{\Upsilon }^L_C$ are the forward generalization of the embedding potentials $\tilde{Y}^{L,x}_C$. Note that, in the model of Meir and Wingreen, the leads are non-interacting; hence, the second trace Tr_λ[...] in equation (17) vanishes.

First we show that the terms in Σ^a/r_L can be factorized out and included within the renormalization of the left lead Green functions g^{a/r, <}_L. Hence, the matrix elements G^<_nλ = 〈n|(GΣg)^<|λ〉 can be recast as $G^<_{n\lambda }=\langle n\vert (G_C\ \Sigma _{CL}\ \tilde{g}_L)^<\vert \lambda \rangle$ with $\tilde{g}_L$ the renormalized L lead Green function whose definition is given below.

We start from

$$\begin{eqnarray} &&\fl G^<_{n\lambda }=\langle n\vert (G\Sigma g)^<\vert \lambda \rangle \nonumber\\ &&=\bullet \hspace{-15.07164pt}\sum _{\lambda _1,\lambda _2,m} G^r_{n\lambda _1}\ \Sigma ^<_{\lambda _1 \lambda _2}\ g^a_{\lambda _2 \lambda } + G^r_{nm}\ \Sigma ^<_{m \lambda _2}\ g^a_{\lambda _2 \lambda } +\ G^<_{n\lambda _1}\ \Sigma ^a_{\lambda _1 \lambda _2}\ g^a_{\lambda _2 \lambda } + G^<_{nm}\ \Sigma ^a_{m \lambda _2}\ g^a_{\lambda _2 \lambda } \nonumber\\ &&+\, G^r_{n\lambda _1}\ \Sigma ^r_{\lambda _1 \lambda _2}\ g^<_{\lambda _2 \lambda } + G^r_{nm}\ \Sigma ^r_{m \lambda _2}\ g^<_{\lambda _2 \lambda } ; \end{eqnarray} \tag{ 30 }$$

hence (we now use Einstein's notation for the sums),

$$\begin{eqnarray} &&\fl G^<_{n\lambda } - G^<_{n\lambda _1}\ \Sigma ^a_{\lambda _1 \lambda _2}\ g^a_{\lambda _2 \lambda } = G^r_{n\lambda _1}\ \Sigma ^<_{\lambda _1 \lambda _2}\ g^a_{\lambda _2 \lambda } + G^r_{nm}\ \Sigma ^<_{m \lambda _2}\ g^a_{\lambda _2 \lambda } + G^<_{nm}\ \Sigma ^a_{m \lambda _2}\ g^a_{\lambda _2 \lambda } \nonumber\\ &&+\, G^r_{n\lambda _1}\ \Sigma ^r_{\lambda _1 \lambda _2}\ g^<_{\lambda _2 \lambda } + G^r_{nm}\ \Sigma ^r_{m \lambda _2}\ g^<_{\lambda _2 \lambda } \end{eqnarray} \tag{ 31 }$$

and

$$\begin{eqnarray} &&\fl G^<_{n\lambda _1} (1 - \Sigma ^a\ g^a)_{\lambda _1 \lambda } = G^r_{n\lambda _1}\ \Sigma ^<_{\lambda _1 \lambda _2}\ g^a_{\lambda _2 \lambda } + G^r_{nm}\ \Sigma ^<_{m \lambda _2}\ g^a_{\lambda _2 \lambda } \nonumber\\ &&+\, G^<_{nm}\ \Sigma ^a_{m \lambda _2}\ g^a_{\lambda _2 \lambda } + G^r_{n\lambda _1}\ \Sigma ^r_{\lambda _1 \lambda _2}\ g^<_{\lambda _2 \lambda } + G^r_{nm}\ \Sigma ^r_{m \lambda _2}\ g^<_{\lambda _2 \lambda } , \end{eqnarray} \tag{ 32 }$$

so

$$\begin{eqnarray} &&\fl G^<_{n\lambda } = G^r_{n\lambda _1}\ \Sigma ^<_{\lambda _1 \lambda _2}\ \tilde{g}^a_{\lambda _2 \lambda } + G^r_{nm}\ \Sigma ^<_{m \lambda _2}\ \tilde{g}^a_{\lambda _2 \lambda } + G^<_{nm}\ \Sigma ^a_{m \lambda _2}\ \tilde{g}^a_{\lambda _2 \lambda } \nonumber\\ &&+\, G^r_{n\lambda _1}\ \Sigma ^r_{\lambda _1 \lambda _2}\ g^<_{\lambda _2 \lambda _1} (1 - \Sigma ^a\ g^a)^{-1}_{\lambda _1 \lambda } + G^r_{nm}\ \Sigma ^r_{m \lambda _2}\ g^<_{\lambda _2 \lambda _1} (1 - \Sigma ^a\ g^a)^{-1}_{\lambda _1 \lambda } , \end{eqnarray} \tag{ 33 }$$

where we define the renormalized Green functions $\tilde{g}^a_L$ for the left L electrode as

$$\begin{equation} g^a_{\lambda \lambda ^{\prime }} ( 1- \Sigma ^{{\rm MB},a} g^a )^{-1}_{\lambda ^{\prime } \lambda _1} = \tilde{g}^a_{\lambda \lambda _1} . \end{equation} \tag{ 34 }$$

To solve for equation (33), we also need to know the matrix element $G^r_{n\lambda _1}$:

$$\begin{equation} \fl G^r_{n\lambda _1}=\langle n\vert (G\Sigma g)^r\vert \lambda _1\rangle =\langle n\vert G^r \Sigma ^r g^r \vert \lambda _1\rangle = G^r_{n\lambda _2}\ \Sigma ^r_{\lambda _2 \lambda _3}\ {g}^r_{\lambda _3 \lambda _1} + G^r_{nm}\ \Sigma ^r_{m \lambda _3}\ {g}^r_{\lambda _3 \lambda _1} ; \end{equation} \tag{ 35 }$$

hence,

$$\begin{equation} G^r_{n\lambda _3} \big(1 -\Sigma ^r_L\ {g}^r_L\big)_{\lambda _3 \lambda _1} = G^r_{nm}\ \Sigma ^r_{m \lambda _3}\ {g}^r_{\lambda _3 \lambda _1} , \end{equation} \tag{ 36 }$$

so

$$\begin{equation} G^r_{n\lambda _1} = G^r_{nm}\ \Sigma ^r_{m \lambda _2}\ \tilde{g}^r_{\lambda _2 \lambda _1} , \end{equation} \tag{ 37 }$$

with a similar definition for the renormalized Green function $\tilde{g}^r_L$ as given for $\tilde{g}^a_L$ in equation (34). With our compact notation, we have

$$\begin{equation} \tilde{g}^{r/a}_L = {g}^{r/a}_L + {g}^{r/a}_L\ \Sigma ^{{\rm MB},r/a}_L\ \tilde{g}^{r/a}_L . \end{equation} \tag{ 38 }$$

Using the result of equation (37) into equation (33) and using the fact that $\big(1 - \Sigma ^a_L\ g^a_L\big)^{-1} = \big(1 + \Sigma ^a_L\ \tilde{g}^a_L\big)$, we find that

$$\begin{eqnarray} &&\fl G^<_{n\lambda } = G^r_{nm}\ \Sigma ^r_{m \lambda _2}\ \tilde{g}^r_{\lambda _2 \lambda _1}\ \Sigma ^<_{\lambda _1 \lambda _2}\ \tilde{g}^a_{\lambda _2 \lambda } + G^r_{nm}\ \Sigma ^<_{m \lambda _2}\ \tilde{g}^a_{\lambda _2 \lambda } + G^<_{nm}\ \Sigma ^a_{m \lambda _2}\ \tilde{g}^a_{\lambda _2 \lambda } \nonumber\\ &&+\, G^r_{nm}\ \Sigma ^r_{m \lambda _2}\ \tilde{g}^r_{\lambda _2 \lambda _1}\ \Sigma ^r_{\lambda _1 \lambda _2}\ g^<_{\lambda _2 \lambda _1} (1 + \Sigma ^a\ \tilde{g}^a)_{\lambda _1 \lambda } + G^r_{nm}\ \Sigma ^r_{m \lambda _2}\ g^<_{\lambda _2 \lambda _1} (1 + \Sigma ^a\ g^a)_{\lambda _1 \lambda } \nonumber\\ &&= G^r_{nm}\ \Sigma ^r_{m \lambda _2}\ \tilde{g}^<_{\lambda _2 \lambda } + G^r_{nm}\ \Sigma ^<_{m \lambda _2}\ \tilde{g}^a_{\lambda _2 \lambda } + G^<_{nm}\ \Sigma ^a_{m \lambda _2}\ \tilde{g}^a_{\lambda _2 \lambda } \nonumber\\ &&= \langle n\vert (G_C\ \Sigma _{CL}\ \tilde{g}_L)^<\vert \lambda \rangle , \end{eqnarray} \tag{ 39 }$$

with

$$\begin{equation} \tilde{g}^<_L = \tilde{g}^{r}_L\ \Sigma ^{{\rm MB},<}_L\ \tilde{g}^{a}_L + \big( 1 + \tilde{g}^{r}_L\ \Sigma ^{{\rm MB},r}_L \big) {g}^<_L \big( 1 + \Sigma ^{{\rm MB},a}_L\ \tilde{g}^{a}_L\big). \end{equation} \tag{ 40 }$$

In other words: the $\tilde{g}^x_L$ are the Green functions in the left L electrode renormalized by the many-body self-energy Σ^{MB, x}_L defined in the same subspace of the left L electrode. Similar results can be derived for the Green functions of the right R electrode.

We find for G^r_C = 〈n|G^r|m〉

$$\begin{eqnarray} G^r_C &=& {g}^r_C + {g}^r_C\ \Sigma ^{{\rm MB},r}_C\ G^r_C + {g}^r_C\ \tilde{Y}^{L+R,r}_C\ G^r_C \nonumber\\ &=& \big[ \big({g}^r_C\big)^{-1} - \Sigma ^{{\rm MB},r}_C - \tilde{Y}^{L+R,r}_C \big]^{-1} \nonumber\\ &=& \big[ \big(\tilde{g}^r_C\big)^{-1} - \tilde{Y}^{L+R,r}_C \big]^{-1} , \end{eqnarray} \tag{ 41 }$$

where $\tilde{Y}^{L+R,r}_C$ is the sum $\tilde{Y}^{L+R,r}_C=\tilde{Y}^{L,r}_C+\tilde{Y}^{R,r}_C$ of the generalized lead self-energies $\tilde{Y}^{\alpha ,r}_C$ (α = L, R) defined previously as $\tilde{Y}^{\alpha ,r}_C = (\Sigma _{C\alpha } \tilde{g}_{\alpha } \Sigma _{\alpha C})^r = \Sigma ^r_{C\alpha } \tilde{g}^r_{\alpha }\Sigma ^r_{\alpha C}$ with Σ_Cα = V_Cα + Σ^MB_Cα.

A similar expression can be derived for the advanced Green function G^a_C in the central region by swapping r↔a.

We find for G^r_C = 〈n|G^<|m〉 that

$$\begin{equation} G^<_C = G^r_C\ \big( \Sigma ^{{\rm MB},<}_C + \tilde{Y}^{L+R,<}_C \big) G^a_C, \end{equation} \tag{ 42 }$$

with $\tilde{Y}^{L+R,<}_C(\omega ) = \tilde{Y}^{L,<}_C(\omega ) + \tilde{Y}^{R,<}_C(\omega ) = \bullet \hspace{-8.53581pt}\sum _{\alpha =L,R} \left( \Sigma _{C\alpha }(\omega ) \tilde{g}_{\alpha }(\omega ) \Sigma _{\alpha C}(\omega ) \right)^<$. Using the rules of analytical continuation for products given in appendix B, one gets $\tilde{Y}^{L+R,<}_C=\bullet \hspace{-8.53581pt}\sum _{\alpha =L,R} \Sigma ^r_{C\alpha } \tilde{g}^r_{\alpha } \Sigma ^<_{\alpha C} + \Sigma ^r_{C\alpha } \tilde{g}^<_{\alpha } \Sigma ^a_{\alpha C} + \Sigma ^<_{C\alpha } \tilde{g}^a_{\alpha } \Sigma ^a_{\alpha C}$.

We find for G^r_R = 〈ρ|G^r|ρ'〉 that

$$\begin{equation} G^r_{\rho \rho ^{\prime }} = \tilde{g}^r_{\rho \rho ^{\prime }} + \tilde{g}^r_{\rho \rho _1}\ \tilde{Y}^{CL,r}_{\rho _1\rho _2} G^r_{\rho _2 \rho ^{\prime }}, \end{equation} \tag{ 43 }$$

where $\tilde{ Y}^{CL,r}_R$ is an embedding potential of the effects of the central region C (connected to the left lead) on the right lead. It is defined as

$$\begin{equation} \tilde{Y}^{CL,r}_{\rho _1\rho _2}(\omega ) = \Sigma ^r_{\rho _1 m}(\omega )\ \big[ \big[\tilde{g}^r_C(\omega )\big]^{-1} - \tilde{Y}^{L,r}_C(\omega ) \big]^{-1}_{ml} \Sigma ^r_{l \rho _2}(\omega ) , \end{equation} \tag{ 44 }$$

with Σ^r_RC = V_RC + Σ^{MB, r}_RC and similarly for Σ^r_CR. And $\big[ \big[\tilde{g}^r_C(\omega )\big]^{-1} - \tilde{Y}^{L,r}_C(\omega ) \big]^{-1}$ is a retarded Green function of the central region renormalized by the interaction inside C and by the lead self-energy/embedding potential $\tilde{Y}^{L,r}_C$ of the left lead only, with $\tilde{g}^r_C$ defined in equation (41) as $\big(\tilde{g}^r_C\big)^{-1} = \big({g}^r_C\big)^{-1} - \Sigma ^{{\rm MB},r}_C$.

Similarly we can find the expressions for the L lead Green function

$$\begin{eqnarray} &&G^r_L = \tilde{g}^r_L + \tilde{g}^r_L \tilde{Y}^{CR,r}_L G^r_L , \\ &&\tilde{Y}^{CR,r}_L = \Sigma ^r_{LC} \big[ \big(\tilde{g}^r_C\big)^{-1} - \tilde{Y}^{R,r}_C \big]^{-1} \Sigma ^r_{CL} . \end{eqnarray} \tag{ 45 }$$

Finally all the expressions given in this section hold for the advanced G^a_{L, R} by swapping r↔a.

We first concentrate on G^<_L = (GΣg + g)^<_LL. The calculations are rather lengthy, but somehow trivial, and include the derivation of many intermediate Green's functions and many re-factorizations. We find, in agreement with the results of the previous section, that

$$\begin{equation} G^<_L = \big(1 + G^r_L \tilde{Y}^{CR,r}_L \big) \tilde{g}^<_L \big(1 + \tilde{Y}^{CR,a}_L G^a_L\big) + G^r_L \tilde{Y}^{CR,<}_L G^a_L \end{equation} \tag{ 46 }$$

and

$$\begin{eqnarray} &&\tilde{Y}^{CR,x}_L = \big(\Sigma _{LC} \tilde{g}^R_C \Sigma _{CL} \big)^x , \nonumber\\ &&\tilde{g}^{R,r/a}_C = \big[ \big(\tilde{g}^{r/a}_C\big)^{-1} - \tilde{Y}^{R,r/a}_C \big]^{-1} , \\ &&\tilde{g}^{R,<}_C = \big(1 + \tilde{g}^{R,r/a}_C \tilde{Y}^{R,r}_C \big) \tilde{g}^<_C \big(1 + \tilde{Y}^{R,a}_C \tilde{g}^{R,a}_C\big) + \tilde{g}^{R,r}_C \tilde{Y}^{R,<}_C \tilde{g}^{R,a}_C . \nonumber \end{eqnarray} \tag{ 47 }$$

The expression for G^<_R can be obtained from equation (46) by swapping the index L to R and the self-energy $\tilde{Y}^{CR,x}_L$ by $\tilde{Y}^{CL,x}_R$ given in equation (44).

Within the partitioned scheme devised by Meir and Wingreen, i.e. interaction only in the central region C, it can be shown that the trace

$$\begin{equation} {\rm Tr}_n [ \Sigma ^<(\omega ) G^>(\omega ) - \Sigma ^>(\omega ) G^<(\omega ) ] = 0 \end{equation} \tag{ 48 }$$

vanishes for each ω, where Σ(ω) is the total self-energy of the region C: Σ = Σ^MB_C + Σ^L_C + Σ^R_C. The α-lead's self-energy is Σ^α_C = V_Cαg_αV_αC.

Equation (48) is derived from the definition G^{>, <}(ω) = G^r(ω)Σ^{>, <}(ω)G^a(ω) in the region C, and hence (G^r)⁻¹(G^> − G^<)(G^a)⁻¹ = Σ^> − Σ^< = Σ^r − Σ^a = (G^a)⁻¹ − (G^r)⁻¹.

With equation (48), we can derive a condition that must be fulfilled by the interaction self-energy [48] in order to satisfy the current conservation condition I_L + I_R = 0 [48]. This condition is given by

$$\begin{equation} \int {\rm d}\omega \ {\rm Tr}_n \left[ \Sigma ^{{\rm MB}<}_C\ G_C^> - \Sigma ^{{\rm MB}>}_C\ G^<_C \right] = 0 , \end{equation} \tag{ 49 }$$

which means that the integrated collision term must vanish. This is a condition familiarly obtained from a Boltzmann-like treatment of scattering theory [48]. When the interaction self-energy Σ^MB is derived from the so-called Φ-derivable approximation [49, 50, 39, 40], it automatically satisfies the condition given by equation (49).

Equation (49) can also be used as a measure of the accuracy of numerical schemes used to calculate approximately the interaction self-energy . It can also be used as a general constraint equation in the determination a new functional forms for the interaction self-energy.

Now, when the interaction exists throughout the entire system, the derivation described above no longer holds and needs to be generalized to the presence of interaction within the leads and crossing at the left and right contacts.

First we consider the general definition of the lesser and greater Green functions :

$$\begin{equation} G^< = ( 1 + G^r \Sigma ^r ) g^< ( 1 + \Sigma ^a\ G^a ) + G^r\ \Sigma ^<\ G^a . \end{equation} \tag{ 50 }$$

The first term represents the initial conditions g^< before the interaction and the coupling between the different regions are applied.

For the central region, we have chosen the initial condition such as 〈n|g^<|m〉 ≡ 0 (see equation (42)). We could have chosen another initial condition. Such choices have no effects on the steady-state regime when a steady current flows through the central region; however, the initial conditions play an important role in the transient behaviour of the current [51–54].

For the definition of the lesser left- and right-lead Green functions, however, it is not possible to neglect the initial conditions (before full interactions and coupling to the region central are taken into account). It would not be physically correct to ignore the presence of the left and right Fermi seas since they are the thermodynamical limit of the two semi-infinite leads, and act as an electron emitter and collector in our model of a device.

By using the standard Dyson equations for G^a/r, we can recast equation (50) as

$$\begin{equation} G^< = G^r\ ( (g^r)^{-1} g^< (g^a)^{-1} + \Sigma ^{<} )\ G^a = G^r\ \bar{\Sigma }^{<}\ G^a , \end{equation} \tag{ 51 }$$

with $\bar{\Sigma }^{<} = \Sigma ^{<} + \gamma ^<$ and γ^< = (g^r)⁻¹g^<(g^a)⁻¹, and similarly for G^>. Hence, γ^< − γ^> = (g^a)⁻¹ − (g^r)⁻¹ and $\bar{\Sigma }^< - \bar{\Sigma }^> = (G^a)^{-1} - (G^r)^{-1}$. From these properties, it can be easily shown that

$$\begin{equation} {\rm Tr}_{\rm all} [ \bar{\Sigma }^< G^> - \bar{\Sigma }^> G^< ] = 0 \end{equation} \tag{ 52 }$$

for each ω. The trace runs over all indices in the system all ≡ {λ, n, ρ} and the interaction Σ is spread over the whole L, C, R regions. This is a generalization of equation (48).

Because the trace runs over all the three subspaces, we can apply the usual cyclic permutation and recast equation (52) as follows:

$$\begin{equation} -{\rm Tr}_{\rm all} [ (\Sigma G)^< - (G \Sigma )^< ] + {\rm Tr}_{\rm all} [ \gamma ^< G^> - \gamma ^> G^< ] = 0 \end{equation} \tag{ 53 }$$

or equivalently

$$\begin{equation} \int {\rm d}\omega \ {\rm Tr}_{\rm all} [ (\Sigma G)^< - (G \Sigma )^< ] + {\rm Tr}_{\rm all} [ \gamma ^> G^< - \gamma ^< G^> ] = 0. \end{equation} \tag{ 54 }$$

Expanding the trace in the first term over each subspace Tr_all[...] = Tr_λ[..] + Tr_n[...] + Tr_ρ[...], one can easily identify the definition of the left I_L and right I_R currents from Tr_λ[...] and Tr_ρ[...], respectively (see equation (9) above).

Hence, the condition of current conservation I_L + I_R = 0 leads to

$$\begin{equation} \int {\rm d}\omega \ {\rm Tr}_n [ (\Sigma G)^< - (G \Sigma )^< ] + {\rm Tr}_{\rm all} [ \gamma ^> G^< - \gamma ^< G^> ] = 0 . \end{equation} \tag{ 55 }$$

After further manipulation of the trace Tr_n[...] using the rules of analytical continuation and the relationship between the different Green functions and self-energies, we find that the current conservation implies that

$$\begin{equation} \fl \int {\rm d}\omega \ {\rm Tr}_n \big[ \big(\Sigma ^{\rm MB}_C + \tilde{Y}^{L+R}_C\big)^> G^<_C - \big(\Sigma ^{\rm MB}_C + \tilde{Y}^{L+R}_C\big)^< G^>_C \big] + {\rm Tr}_{L,R} [ \gamma ^> G^< - \gamma ^< G^> ] = 0 , \end{equation} \tag{ 56 }$$

where in the second trace the sum runs only over the left and right subspaces, because we have chosen the initial condition for the central region C such that γ^{>, <} = 0.

What is really interesting with the first trace Tr_n[...] in equation (56) is that it has exactly the same form as equation (52), but with all quantities (as well as the trace) defined within the central region C only (i.e. $\bar{\Sigma }\equiv \Sigma ^{\rm MB}_C + \tilde{Y}^{L+R}_C$ in the region C). By looking at the definition of G^{>, <}_C given by equation (42), we can establish that, equivalently to equation (52), the trace actually vanishes:

$$\begin{equation} {\rm Tr}_n \big[ \big(\Sigma ^{\rm MB}_C + \tilde{Y}^{L+R}_C\big)^> G^<_C - \big(\Sigma ^{\rm MB}_C + \tilde{Y}^{L+R}_C\big)^< G^>_C \big]=0 \end{equation} \tag{ 57 }$$

for each ω. Therefore, equation (56) reduces to

$$\begin{equation} \int {\rm d}\omega \ {\rm Tr}_{L,R} [ \gamma ^> G^< - \gamma ^< G^> ] = 0 , \end{equation} \tag{ 58 }$$

a condition which however is almost systematically verified (see appendix C).

Hence, it is better to consider equation (57) to find the condition imposed by the current conservation. Indeed by treating each contribution in equation (57) separately and by integrating over the energy, we can introduce the definition of the currents I_L and I_R such as

$$\begin{equation} \fl \int {\rm d}\omega \ {\rm Tr}_n \left[ \Sigma ^{{\rm MB}>}_C G^<_C - \Sigma ^{{\rm MB}<}_C G^>_C \right] -i_L(\omega ) + \Delta i_L(\omega ) -i_R(\omega ) + \Delta i_R(\omega ) = 0 , \end{equation} \tag{ 59 }$$

where I_α = e/h∫i_α(ω)dω and $\Delta i_\alpha (\omega )= \big[ \tilde{Y}^{\alpha ,>}_C G^<_C - \tilde{Y}^{\alpha ,<}_C G^>_C \big] + i_L(\omega )$.

Hence, the condition of current conservation I_L + I_R = 0 leads to the final important result of this section:

$$\begin{eqnarray} &&\fl \int {\rm d}\omega {\rm Tr}_{\alpha =L,C,R} \left[ \Sigma ^{{\rm MB}>}_\alpha G^<_\alpha - \Sigma ^{{\rm MB}<}_\alpha G^>_\alpha \right] + {\rm Tr}_{n} \big[ \big(\tilde{Y}^{L,>}_C + \tilde{\Upsilon }^L_C - \big(\tilde{\Upsilon }^L_R\big)^\dag - \tilde{\Upsilon }^{L,l}_C\big) G^<_C\nonumber\\ &&-\, \big(\tilde{Y}^{L,<}_C - \tilde{\Upsilon }^{L,l}_C\big) G^>_C + \big(\tilde{\Upsilon }^{L,l}_C + \big(\tilde{\Upsilon }^{L,l}_C\big)^\dag \big) G^a \big] \nonumber\\ &&+\, {\rm Tr}_{n} \left[ \lbrace L \leftrightarrow R \rbrace \right] \nonumber\\ &&= 0 . \end{eqnarray} \tag{ 60 }$$

Using the properties of the $\tilde{\Upsilon }^L_C$ quantities (see section 3.3), we can rewrite the condition imposed by the current conservation as

$$\begin{eqnarray} &&\fl \int {\rm d}\omega \ {\rm Tr}_{\alpha =L,C,R} \left[ \Sigma ^{{\rm MB}>}_\alpha G^<_\alpha - \Sigma ^{{\rm MB}<}_\alpha G^>_\alpha \right] \nonumber\\ &&+\, {\rm Tr}_{n} \big[ \big(\tilde{Y}^{L,>}_C - \tilde{\Upsilon }^{L,g}_C\big) G^<_C - \big(\tilde{Y}^{L,<}_C - \tilde{\Upsilon }^{L,l}_C\big) G^>_C + \big(\tilde{\Upsilon }^{L,l}_C + \big(\tilde{\Upsilon }^{L,l}_C\big)^\dag \big) G^a \big] \nonumber\\ &&+\, {\rm Tr}_{n} \left[ \lbrace L \leftrightarrow R \rbrace \right] \nonumber\\ &&= 0 . \end{eqnarray} \tag{ 61 }$$

To understand fully the conditions of current conservation, we make the following observations.

• (i)  
  Equation (61) is the generalization of equation (49) for the systems where the interaction spreads throughout. Like equation (49), it contains similar terms involving the left L and right R lead as well. But it also contains terms arising from the fact that the interaction is crossing at the LC and RC contacts. However, the physical interpretation of equation (61) still corresponds to the fact that the total integrated collision terms must vanish.
• (ii)  
  For interaction present only within the C region, one can show that $\tilde{\Upsilon }^{L,l}_C + \big(\tilde{\Upsilon }^{L,l}_C\big)^\dag =0$ as well as $\tilde{Y}^{L,<}_C - \tilde{\Upsilon }^{L,l}_C = 0$ and $\tilde{Y}^{L,>}_C - \tilde{\Upsilon }^{L,g}_C = 0$ since Σ^MB_LC = V_LC, and likewise for the terms involving the RC interface. Furthermore, in that case, there are no interactions in the leads Σ^MB_{α = L, R} = 0, and hence one recovers equation (49) as expected.
• (iii)  
  Equation (61) is also consistent with the one of the main point made in section 3.3: the current I_α (α = L or R) cannot be obtained from a trace over the central region defined as Tr_n[Y^{α, <}_CG^>_C − Y^{α, >}_CG^<_C]. If it were the case, then the Δi_α(ω) defined after equation (59) is $\Delta i_\alpha (\omega )= \left[ \tilde{Y}^{\alpha ,>}_C G^<_C - \tilde{Y}^{\alpha ,<}_C G^>_C \right] + i_L(\omega ) \equiv 0$. And the current conservation condition would imply that equation (61) reduces to

  $$\begin{equation} \int {\rm d}\omega \ {\rm Tr}_{\alpha =\lambda ,n,\rho } \left[ \Sigma ^{{\rm MB}>}_\alpha G^<_\alpha - \Sigma ^{{\rm MB}<}_\alpha G^>_\alpha \right] = 0. \end{equation} \tag{ 62 }$$

  This equation concerns the collision terms within each of the three regions, but is not complete because it does not show any constraint on the interaction crossing at the LC and RC contacts.

There is one arbitrary choice in our derivation: the location of the LC and RC interfaces with respect to the physical realistic system. Such locations are somehow arbitrary in our formalism but could be conveniently chosen for practical numerical calculations. We have already seen that, in some cases, a local quasi-equilibrium is reached within the left and right leads. In such cases, we get simplified results for the current, since the local non-equilibrium distribution functions f^<_{L, R}, f^{int <}_{L, R} are equal to the local Fermi distribution in the left or right lead f_{L, R} = f^{0 <}_{L, R}.

Our formalism provides a formal justification of the concept of the extended molecule that has been used so far within the conventional partitioned scheme with interaction present only in the C region. Being general, our formalism also provides the corresponding 'correction' terms needed when the interactions cross at the contacts and when the contacts are not in their respective local (quasi) equilibrium.

Our work also justifies the scheme recently used in [45] where two closed surfaces are used to define two sets of LC and RC interfaces around the molecule. Within the inner region, the electron–electron interaction is calculated via the GW scheme with the many-body self-energy Σ(τ, τ') = G(τ, τ')W(τ, τ'). While the screened Coulomb interaction W = v + vPW and the polarization P are calculated for the extended region in order to ensure a better treatment of non-local screening effects. To some extent, this is an empirical way of defining the concept of generalized embedding potentials that we obtain in a formal manner in our formalism.

An extension of the present formalism to non-orthogonal basis sets as often used in electronic structure calculations with localized basis sets would potentially be useful. This has already been achieved for interaction present only in the central region by using a representation in the dual basis in which the lead subspaces are orthogonal to the central region subspace [30].

In a broader context, our formalism introduces in a formal manner the generalization of the concept embedding potential to interacting cases. In the case originally studied by Meir and Wingreen, where the interaction is present only in the central region, the effects of the leads on the Green functions defined with the subspace of the region C are taken into account via the so-called lead self-energies : Σ^α_C(ω) = V_Cαg_α(ω)V_αC. These non-local self-energies are simply a matrix representation of the so-called embedding potential originally defined in real space by Inglesfield [41, 76, 77, 42, 43]. The latter can be seen as defining a surface (or two interfaces) around the central region that the interaction does not cross. The non-locality of the embedding potentials arises only from the Green functions defined on this surface.

In our formalism, when the interaction can cross the LC/RC interfaces, we obtain in a systematic way a generalization of the embedding potentials $\tilde{Y}^\alpha _C$, defined as $\tilde{Y}^\alpha _C(\omega ) = \Sigma _{C\alpha }(\omega ) \tilde{g}_{\alpha }(\omega ) \Sigma _{\alpha C}(\omega )$. These generalized embedding potentials contain a 'double' non-locality, in the sense that the Σ^MB_αC part of Σ_αC can have a larger spatial extent than the hopping matrix elements V_αC. Hence the $\tilde{Y}^\alpha _C$ defines somehow not a simple surface but a 'buffer' zone that is contained between two surfaces whose separation is related to the characteristic (spatial decay) length of the interaction self-energies $\Sigma ^{\rm MB}_{\alpha C} \equiv \Sigma ^{\rm MB}(\vert \mathbf {x}_\alpha -\mathbf {x}_n \vert )$.

For practical numerical calculations of the current, one needs to choose the form for the interaction self-energies. They can be obtained from conventional many-body perturbation theory and Feynman diagrammatics, extended onto the Keldysh contour. The interaction self-energy Σ^MB can be obtained from the Φ-derivable conserving approximation [49, 50, 39, 40], and then they should automatically satisfy the condition of current conservation. We have given an example of such interaction self-energies in [35] for the case of electron–phonon coupling inside the region C and crossing at the LC interface. One could also devise other functional forms for the self-energies, such as the functional of the charge and spin densities, and or of the current density itself. In these cases, one should devise the functionals such that the condition of current conservation is indeed fulfilled.

Finally, one should note that once the functional forms of the self-energies are chosen, one can perform the calculations self-consistently. The self-energies in the three regions L, C, R and at the interfaces LC and RC are functionals of the other electron (and phonon) Green functions (or other physical quantities related to them such as the charge, spin or current density) defined inside the L, C, R regions as well as at the two LC and RC interfaces. Hence the Green functions and self-energies need to be determined self-consistently in all the parts of the system in order to get the current of a fully many-body interacting nanojunction at non-equilibrium.

When modelling the self-energies such that $\tilde{\Upsilon }^{L,l}_C+\big(\tilde{\Upsilon }^{L,l}_C\big)^\dag =0$, we can express the current equation (17) as

$$\begin{eqnarray} \fl I_L = \frac{e}{\hbar } \int \frac{{\rm d}\omega }{2\pi } & {\rm Tr}_n\big[ \big(G^r_C-G^a_C\big) \tilde{\Upsilon }^{L,l}_C + G^<_C\big(\tilde{\Upsilon }^L_C - (\tilde{\Upsilon }^L_C)^\dagger \big) \big] + {\rm Tr}_\lambda \left[\Sigma ^{{\rm MB},>}_{L} G^<_{L} - \Sigma ^{{\rm MB},<}_{L} G^>_{L} \right] , \end{eqnarray} \tag{ 63 }$$

which bear even more resemblance to the Meir and Wingreen result, equation (28), and hence could be recast as a generalized Landauer-like expression for the current [14]. In such cases, the condition of current conservation becomes

$$\begin{eqnarray} &&\fl \int {\rm d}\omega \ {\rm Tr}_{\alpha =L,C,R} \left[ \Sigma ^{{\rm MB}>}_\alpha G^<_\alpha - \Sigma ^{{\rm MB}<}_\alpha G^>_\alpha \right] + {\rm Tr}_{n} \big[ \big(\tilde{Y}^{L,>}_C - \tilde{\Upsilon }^{L,g}_C\big) G^<_C - \big(\tilde{Y}^{L,<}_C - \tilde{\Upsilon }^{L,l}_C\big) G^>_C \big] \nonumber\\ &&+\, {\rm Tr}_{n} \left[ \lbrace L \leftrightarrow R \rbrace \right] \nonumber\\ &&= 0 . \end{eqnarray} \tag{ 64 }$$

Equation (64) looks just like the sum of the integrated collision terms in the form of

$$\begin{equation} \int {\rm d}\omega \ {\rm Tr}_{\beta } \left[ \Sigma ^{{\rm int} >}_\beta G^<_\beta - \Sigma ^{{\rm int} <}_\beta G^>_\beta \right] = 0, \end{equation} \tag{ 65 }$$

where Σ^int_β represents some interaction self-energy and the sum β is over not only the three subspaces L, C and R but also the two interfaces LC and RC.

By using the rules of analytical continuation and the relationships between the different Green functions , one can find that, in general, the sum $\tilde{\Upsilon }^l_{LC}+(\tilde{\Upsilon }^l_{LC})^\dag$ is given by

$$\begin{eqnarray} &&\fl \tilde{\Upsilon }^{L,l}_C + \big(\tilde{\Upsilon }^{L,l}_C\big)^\dag = (\Sigma \tilde{g})^<_{CL}\ \Sigma ^>_{LC} - (\Sigma \tilde{g})^>_{CL}\ \Sigma ^<_{LC} + \Sigma ^>_{CL}\ (\tilde{g} \Sigma )^<_{LC} - \Sigma ^<_{CL}\ (\tilde{g} \Sigma )^>_{LC} . \end{eqnarray} \tag{ 66 }$$

There are two different cases in which the sum $\tilde{\Upsilon }^{L,l}_C + \big(\tilde{\Upsilon }^{L,l}_C\big)^\dag$ vanishes:

• (i)  
  when the interaction is only present in the central region. We have already discussed that case in length in the paper.;
• (ii)  
  when interaction is instantaneous, i.e. local in time Σ^MB(τ, τ') = v^MB(τ)δ(τ − τ'). Hence, there are no lesser/greater components of the self-energy Σ^{MB, > <} = 0. This occurs, as already discussed, in mean-field-based and in density-functional-based theories for which a single (quasi-)particle description of the system is available.

Finally, it would be interesting to study and find cases which go beyond mean-field-based or density-functional-based methods, and for which $\tilde{\Upsilon }^{L,l}_C + \big(\tilde{\Upsilon }^{L,l}_C\big)^\dag =0$. If such cases exist, it would still be possible to analyse their transport properties in terms of a generalized Landauer-like approach [14].