Nonequilibrium perturbation theory in Liouville–Fock space for inelastic electron transport

We begin by considering a quantum system (e.g. quantum dot, molecule, etc) connected to two electrodes, left and right, with different chemical potentials. Each electrode is partitioned into two parts (figure 1): the macroscopically large lead (environment) and the finite buffer zone between the system and the environment. So the Hamiltonian of the whole system is

$$\begin{eqnarray} &&\mathcal{H}={H}_{\mathrm{S}}+{H}_{\mathrm{SB}}+{H}_{\mathrm{B}}+{H}_{\mathrm{BE}}+{H}_{\mathrm{E}}. \end{eqnarray} \tag{ 1 }$$

We assume that the environment and the buffer zones are described by the noninteracting Hamiltonians

$$\begin{eqnarray} &&{H}_{\mathrm{E}}=\mathop{\sum }\limits _{k \alpha }{\varepsilon }_{k \alpha }a_{k \alpha }^{\dagger }{a}_{k \alpha },\tqs {H}_{\mathrm{B}}=\mathop{\sum }\limits _{b \alpha }{\varepsilon }_{b \alpha }a_{b \alpha }^{\dagger }{a}_{b \alpha }. \end{eqnarray} \tag{ 2 }$$

Here ε_kα denote the continuum single-particle spectra of the left (α = L) and right (α = R) lead states, and ${a}_{k \alpha }^{\dagger }~({a}_{k \alpha })$ create (annihilate) an electron in the lead state kα. The buffer zones have discrete energy spectrum ε_bα with corresponding creation and annihilation operators ${a}_{b \alpha }^{\dagger }$ and a_bα. The system Hamiltonian is taken in the most general form:

$$\begin{eqnarray} &&{H}_{\mathrm{S}}=\mathop{\sum }\limits _{s}{\varepsilon }_{s}a_{s}^{\dagger }{a}_{s}+{H}_{\mathrm{S}}^{\prime}, \end{eqnarray} \tag{ 3 }$$

where ${a}_{s}^{\dagger }~({a}_{s})$ create (annihilate) electrons in the single-particle state ε_s in the dot and ${H}_{\mathrm{S}}^{\prime}$ contains two-particle electron–electron correlations and/or electron–vibration coupling. The buffer–environment and quantum dot–buffer couplings have the standard tunneling form:

$$\begin{eqnarray} {H}_{\mathrm{BE}}&=&\mathop{\sum }\limits _{b k \alpha }({v}_{b k \alpha }a_{b \alpha }^{\dagger }{a}_{k \alpha }+\mathrm{h.c.}), \end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray} {H}_{\mathrm{SB}}&=&\mathop{\sum }\limits _{s b \alpha }({t}_{s b \alpha }a_{b \alpha }^{\dagger }{a}_{s}+\mathrm{h.c.}). \end{eqnarray} \tag{ 5 }$$

Zoom In Zoom Out Reset image size

Now we introduce an embedded system which consists of the quantum system itself and the buffer zones. We have recently demonstrated that if we take the buffer zones sufficiently large the density matrix of the embedded system obeys the kinetic equation of Lindblad type. The technical details of the derivations and underlying approximations can be found in the appendix to [24]. Here we give only a sketch of the derivation with the emphasis on important physics relevant to our subsequent discussion.

The starting point is the Liouville equation for the total density matrix χ(t) in the interaction picture:

$$\begin{eqnarray} &&{\dot {\chi }}_{I}(t)=-\mathrm{i}[{v}_{I}(t),{\chi }_{I}(t)]. \end{eqnarray} \tag{ 6 }$$

Here the buffer–environment coupling H_BE is treated as an interaction Hamiltonian, i.e.  = h + H_BE and v_I(t) = e^ihtH_BEe^−iht. To derive the Lindblad master for the reduced density matrix of the embedded system, ρ_I(t) = Tr_Eχ_I(t), we take the trace over the environment in equation (6) and make the following approximations.

• (i)  
  The total density matrix can be factorized as χ_I(t) = ρ_I(t)ρ_E, where ρ_E is the density matrix of the environment taken in the equilibrium grand canonical ensemble form (Born approximation).
• (ii)  
  The environment relaxation time is very fast, so we can use the local-time (Markov) approximation for the reduced density matrix.
• (iii)  
  The single-particle states in the buffer zone propagate as free states:

  $$\begin{eqnarray} &&{\mathrm{e}}^{\mathrm{i}h t}{a}_{b \alpha }{\mathrm{e}}^{-\mathrm{i}h t}={\mathrm{e}}^{-\mathrm{i}{\varepsilon }_{b \alpha }t}{a}_{b \alpha }+\mathrm{O}(1/\sqrt{N}) \end{eqnarray} \tag{ 7 }$$

  where N is the number of discrete single-particle levels of the buffer zone.
• (iv)  
  Rapidly oscillating terms proportional to exp[i(ε_bα − ε_b'α)] for ε_bα ≠ ε_b'α are neglected (rotating wave approximation).

Under these approximations, the Liouville equation (6) reduces to a master equation for the reduced density matrix in Lindblad form. In the Schrödinger representation it can be written as

$$\begin{eqnarray} &&\frac{\mathrm{d}\rho (t)}{\mathrm{d}t}=-\mathrm{i}[H,\rho (t)]+\Pi \rho (t). \end{eqnarray} \tag{ 8 }$$

Here the Hamiltonian H includes the Lamb shift of the single-particle levels of the buffer zones:

$$\begin{eqnarray} &&H={H}_{\mathrm{S}}+{H}_{\mathrm{SB}}+{H}_{\mathrm{B}}+\mathop{\sum }\limits _{b \alpha }{\Delta }_{b \alpha }a_{b \alpha }^{\dagger }{a}_{b \alpha }, \end{eqnarray} \tag{ 9 }$$

and the non-Hermitian dissipator is given by the standard Lindblad form:

$$\begin{eqnarray} &&\fl \Pi \rho (t)=\mathop{\sum }\limits _{b \alpha }\mathop{\sum }\limits _{\mu =1,2}(2{L}_{b \alpha \mu }\rho (t)L_{b \alpha \mu }^{\dagger }-\{ L_{b \alpha \mu }^{\dagger }{L}_{b \alpha \mu },\rho (t)\} ). \end{eqnarray} \tag{ 10 }$$

The operators L_bα1 and L_bα2 are referred to as the Lindblad operators, which represent the buffer–environment interaction. They have the following form:

$$\begin{eqnarray} &&{L}_{b \alpha 1}=\sqrt{{\Gamma }_{b \alpha 1}}{a}_{b \alpha },\tqs {L}_{b \alpha 2}=\sqrt{{\Gamma }_{b \alpha 2}}{a}_{b \alpha }^{\dagger } \end{eqnarray} \tag{ 11 }$$

with Γ_bα1 = γ_bα(1 − f_bα),Γ_bα2 = γ_bαf_bα. Here f_bα = [1+e^{β(ε_bα−μ_α)}]⁻¹ and γ_bα (Δ_bα) is the imaginary (real) part of the environment self-energy ∑_k|v_bkα|²/(ε_bα − ε_kα + i0⁺).

The Lindblad master equation describes the time evolution of the open embedded quantum system preserving the probability and the positivity of the density matrix. Open boundary conditions are taken into account by the non-Hermitian dissipative part of equation (8), Πρ(t), which represents the influence of environment on the buffer zone. The applied bias potential enters into equation (8) via fermionic occupation numbers f_bα which depend on the temperature (β = 1/T) and the chemical potential μ_α in the left and right electrodes.

Let us convert the Lindblad master equation (8) to a non-Hermitian field theory suitable for perturbative many-body calculations. To this aim we need to introduce the concept of creation and annihilation superoperators acting on the Liouville–Fock space [25–27, 22]. Our introduction of the Liouville–Fock space closely follows the work of Schmutz [25]. It is general and not restricted to a particular choice of the kinetic equation.

Let {|n)} be a complete orthonormal basis set in the Fock space :

$$\begin{eqnarray} &&\mathop{\sum }\limits _{n}\vert n)(n\vert =I,\tqs (n\vert m)={\delta }_{n m}. \end{eqnarray} \tag{ 12 }$$

It is formed by particle number eigenstates $\vert n)={a}_{{j}_{1}}^{\dagger }\cdots {a}_{{j}_{n}}^{\dagger }\vert 0)$, such that ${a}_{j}^{\dagger }{a}_{j}\vert n)={n}_{j}\vert n)$. Here |0) is the vacuum state and ${a}_{j}^{\dagger },{a}_{j}$ are creation and annihilation operators for single-particle state j. Without loss of generality we focus on fermions, so we assume that ${a}_{j}^{\dagger }$ and a_j satisfy the canonical anti-commutation relations.

The set of linear operators {A(a^†,a)} acting on form a linear vector space, which is called the Liouville–Fock space associated with . We denote an element of the Liouville–Fock space by  |A〉. The scalar product of two elements of the Liouville–Fock space is defined as

$$\begin{eqnarray} &&\langle {A}_{1}\vert {A}_{2}\rangle =\mathrm{Tr}({A}_{1}^{\dagger }{A}_{2}). \end{eqnarray} \tag{ 13 }$$

In the Liouville–Fock space we introduce a complete orthonormal basis {|m,n〉 =  ||m)(n|〉}, which satisfies

$$\begin{eqnarray} &&\langle m n\vert {m}^{\prime}{n}^{\prime}\rangle ={\delta }_{m{m}^{\prime}}{\delta }_{n{n}^{\prime}},\tqs \mathop{\sum }\limits _{m n}\hspace{0.167em} \vert m n\rangle \langle m n\vert \hspace{0.167em} =\bar {I}. \end{eqnarray} \tag{ 14 }$$

Here 〈mn|  =  |mn〉^† = 〈[|m)(n|]^†|  = 〈|n)(m||  and $\bar {I}$ is the identity operator in the Liouville–Fock space. Then, for an arbitrary element of the Liouville–Fock space we have

$$\begin{eqnarray} &&\hspace{0.167em} \vert A\rangle =\mathop{\sum }\limits _{m n}{A}_{m n}\hspace{0.167em} \vert m n\rangle , \end{eqnarray} \tag{ 15 }$$

where A_mn = 〈m|A|n〉 = 〈mn|A〉. In particular, the identity operator I in equation (12) corresponds to

$$\begin{eqnarray} &&\hspace{0.167em} \vert I\rangle =\mathop{\sum }\limits _{n}\hspace{0.167em} \vert n,n\rangle . \end{eqnarray} \tag{ 16 }$$

The scalar product of a vector  |A〉 with 〈I|  is equivalent to the trace operation in the Fock space:

$$\begin{eqnarray} &&\langle I\vert A\rangle =\mathrm{Tr}(A), \end{eqnarray} \tag{ 17 }$$

and for the density matrix we have 〈I|ρ〉 = 1.

As was suggested by Schmutz [25] we introduce superoperators $\hat {a},\tilde {a}$ through their action on the basis vectors  |mn〉:

$$\begin{eqnarray} &&\fl {\hat {a}}_{j}\hspace{0.167em} \vert m n\rangle =\hspace{0.167em} \vert {a}_{j}\vert m)(n\vert \rangle ,\tqs {\tilde {a}}_{j}\hspace{0.167em} \vert m n\rangle =\mathrm{i}(-1)^{\mu }\hspace{0.167em} \vert \vert m)(n\vert {a}_{j}^{\dagger }\rangle ,\nonumber\\ \end{eqnarray} \tag{ 18 }$$

where μ = ∑_j(m_j + n_j) = m + n. By analyzing the Hermitian conjugate of the matrix elements of $\hat {a},\tilde {a}$, we find

$$\begin{eqnarray} &&\fl {\hat {a}}_{j}^{\dagger }\hspace{0.167em} \vert m n\rangle =\hspace{0.167em} \vert {a}_{j}^{\dagger }\vert m)(n\vert \rangle ,\tqs {\tilde {a}}_{j}^{\dagger }\hspace{0.167em} \vert m n\rangle =\mathrm{i}(-1)^{\mu }\hspace{0.167em} \vert \vert m)(n\vert {a}_{j}\rangle .\nonumber\\ \end{eqnarray} \tag{ 19 }$$

It follows from (18) and (19) that superoperators $\hat {a},{\hat {a}}^{\dagger }$ simulate the action of a and a^† on |m)(n| from the left, while $\tilde {a},{\tilde {a}}^{\dagger }$ simulate the action of a^† and a on |m)(n| from the right. Here we would like to emphasize that our definition of tilde superoperators $\tilde {a},{\tilde {a}}^{\dagger }$ differs from Schmutz's definition by phase factors  − i and  + i, respectively. The reason for introducing these factors is that the so-called tilde substitution rule (see below) becomes simpler. We also note that the alternative definition for superoperators is used in [27], where the 'right' creation and annihilation superoperators are not Hermitian conjugate to each other.

As follows from (18) and (19), the superoperators ${\hat {a}}_{j},{\tilde {a}}_{j},{\hat {a}}_{j}^{\dagger },{\tilde {a}}_{j}^{\dagger }$ obey the fermionic anti-commutation relations:

$$\begin{eqnarray} &&\{ {\hat {a}}_{i},{\hat {a}}_{j}^{\dagger }\} =\{ {\tilde {a}}_{i},{\tilde {a}}_{j}^{\dagger }\} ={\delta }_{i j}, \end{eqnarray} \tag{ 20 }$$

while other anti-commutators vanish:

$$\begin{eqnarray} &&\{ {\hat {a}}_{i},{\hat {a}}_{j}\} =\{ {\tilde {a}}_{i},{\tilde {a}}_{j}\} =\{ {\hat {a}}_{i},{\tilde {a}}_{j}\} =\{ {\hat {a}}_{i},{\tilde {a}}_{j}^{\dagger }\} =0. \end{eqnarray} \tag{ 21 }$$

It also follows from (18) and (19) that $\hat {a}\hspace{0.167em} \vert 0 0\rangle =\tilde {a}\hspace{0.167em} \vert 0 0\rangle =0$ and the Liouville–Fock space basis vectors are generated from the vacuum  |00〉 by application of the creation superoperators

$$\begin{eqnarray} &&\hspace{0.167em} \vert m n\rangle =(-\mathrm{i})^{{n}^{2}}{\hat {a}}_{{k}_{1}}^{\dagger }\cdots {\hat {a}}_{{k}_{m}}^{\dagger }{\tilde {a}}_{{l}_{1}}^{\dagger }\cdots {\tilde {a}}_{{l}_{n}}^{\dagger }\hspace{0.167em} \vert 0 0\rangle . \end{eqnarray} \tag{ 22 }$$

Moreover, basis vectors  |mn〉 are 'super-fermion' number eigenstates:

$$\begin{eqnarray} &&{\hat {a}}_{j}^{\dagger }{\hat {a}}_{j}\hspace{0.167em} \vert m n\rangle ={m}_{j}\hspace{0.167em} \vert m n\rangle ,\tqs {\tilde {a}}_{j}^{\dagger }{\tilde {a}}_{j}\hspace{0.167em} \vert m n\rangle ={n}_{j}\hspace{0.167em} \vert m n\rangle . \end{eqnarray} \tag{ 23 }$$

Using the definition of superoperators we can rewrite the identity (16) in the following form:

$$\begin{eqnarray} &&\hspace{0.167em} \vert I\rangle =\exp \left (-\mathrm{i}\mathop{\sum }\limits _{j}\hat {a}_{j}^{\dagger }{\tilde {a}}_{j}^{\dagger }\right )\hspace{0.167em} \vert 0 0\rangle . \end{eqnarray} \tag{ 24 }$$

Note that, because of the different definition of tilde superoperators, the obtained expression for  |I〉 differs from Schmutz's analogous expression [25] by the phase factor (−i) in the exponent. From (18), (19) and (24) we find that the superoperators ${\hat {a}}_{j}^{\dagger }$ and ${\hat {a}}_{j}$ are connected to their tilde conjugate ${\tilde {a}}_{j}^{\dagger }$ and ${\tilde {a}}_{j}$ by the relations

$$\begin{eqnarray} &&{\hat {a}}_{j}\hspace{0.167em} \vert I\rangle =-\mathrm{i}{\tilde {a}}_{j}^{\dagger }\hspace{0.167em} \vert I\rangle ,\tqs {\hat {a}}_{j}^{\dagger }\hspace{0.167em} \vert I\rangle =-\mathrm{i}{\tilde {a}}_{j}\hspace{0.167em} \vert I\rangle . \end{eqnarray} \tag{ 25 }$$

For an operator A = A(a^†,a) given by the power series of creation and annihilation operators we define two superoperators:

$$\begin{eqnarray} &&\hat {A}=A({\hat {a}}^{\dagger },\hat {a}),\tqs \tilde {A}={A}^{\ast }({\tilde {a}}^{\dagger },\tilde {a}). \end{eqnarray} \tag{ 26 }$$

Here, the ∗ means the complex conjugate of the c-number coefficients. The relation between non-tilde and tilde superoperators is given by the following tilde conjugation rules:

$$\begin{eqnarray} && \displaystyle ({c}_{1}{\hat {A}}_{1}+{c}_{2}{\hat {A}}_{2})^{\widetilde {}}={c}_{1}^{\ast }{\tilde {A}}_{1}+{c}_{2}^{\ast }{\tilde {A}}_{2},\\ && \displaystyle ({\hat {A}}_{1}{\hat {A}}_{2})^{\widetilde {}}={\tilde {A}}_{1}{\tilde {A}}_{2},\tqs (\tilde {A})^{\widetilde {}}=A. \end{eqnarray} \tag{ 27 }$$

Applying tilde conjugation to  |mn〉 we find

$$\begin{eqnarray} &&\hspace{0.167em} \vert m n\rangle ^{\widetilde {}}=(+\mathrm{i})^{{\mu }^{2}}\hspace{0.167em} \vert n m\rangle , \end{eqnarray} \tag{ 28 }$$

where μ = m + n. Therefore $\hspace{0.167em} \vert I\rangle ^{\widetilde {}}=\hspace{0.167em} \vert I\rangle$, i.e.  |I〉 is tilde-invariant. Generally, if A = A(a^†,a) is a Hermitian bosonic operator then $\hspace{0.167em} \vert A\rangle ^{\widetilde {}}=\hspace{0.167em} \vert A\rangle$.

According to the definition of the superoperator $\hat {A}$, if A = ∑_mnA_mn|m)(n| then $\hat {A}={\sum }_{m n k}{A}_{m n}\hspace{0.167em} \vert m k\rangle \langle n k\vert \hspace{0.167em}$ and we obtain

$$\begin{eqnarray} &&\hspace{0.167em} \vert A\rangle =\hat {A}\hspace{0.167em} \vert I\rangle ,\tqs \hspace{0.167em} \vert A\rangle ^{\widetilde {}}=\tilde {A}\hspace{0.167em} \vert I\rangle , \end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} &&\hspace{0.167em} \vert {A}_{1}{A}_{2}\rangle ={\hat {A}}_{1}{\hat {A}}_{2}\hspace{0.167em} \vert I\rangle ={\hat {A}}_{1}\hspace{0.167em} \vert {A}_{2}\rangle . \end{eqnarray} \tag{ 30 }$$

Therefore, the expectation value of an operator A = A(a^†,a) in the state with the density matrix ρ = ρ(a^†,a) can be calculated as the matrix element of the corresponding superoperator $\hat {A}=A({\hat {a}}^{\dagger },\hat {a})$ sandwiched between 〈I|  and $\hspace{0.167em} \vert \rho \rangle =\hat {\rho }\hspace{0.167em} \vert I\rangle$:

$$\begin{eqnarray} &&\langle A\rangle =\mathrm{Tr}(A \rho )=\langle I\vert A \rho \rangle =\langle I\vert \hat {A}\vert \rho \rangle . \end{eqnarray} \tag{ 31 }$$

Using (25) we can show that the following tilde substitution rule is valid:

$$\begin{eqnarray} &&\hat {A}\hspace{0.167em} \vert I\rangle ={\sigma }_{A}{\tilde {A}}^{\dagger }\hspace{0.167em} \vert I\rangle . \end{eqnarray} \tag{ 32 }$$

Here σ_A =+ 1 if A is a bosonic operator and σ_A =− i if A is a fermionic operator. Moreover, taking into account that non-tilde and tilde fermion superoperators anti-commute we find that

$$\begin{eqnarray} &&{\hat {A}}_{1}\hspace{0.167em} \vert {A}_{2}\rangle =\mathrm{i}{\tilde {A}}_{2}^{\dagger }\hspace{0.167em} \vert {A}_{1}\rangle , \end{eqnarray} \tag{ 33 }$$

if both A₁ and A₂ are fermionic operators, and

$$\begin{eqnarray} &&{\hat {A}}_{1}\hspace{0.167em} \vert {A}_{2}\rangle ={\sigma }_{{A}_{2}}{\tilde {A}}_{2}^{\dagger }\hspace{0.167em} \vert {A}_{1}\rangle \end{eqnarray} \tag{ 34 }$$

otherwise. It should be noted that the Schmutz tilde substitution rule [25] is cumbersome and it takes the simple form like (32) only if all terms in the power series of A(a^†,a) have the common quantity m − n. Here m(n) is the number of creation (annihilation) operators.

The general prescription to obtain the equation for  |ρ(t)〉 from the kinetic equation for ρ(t) is the following. First, we transform the kinetic equation for ρ = ρ(a^†,a) into the kinetic equation for $\hat {\rho }=\rho ({\hat {a}}^{\dagger },\hat {a})$ by formally replacing all operators a^†,a by superoperators ${\hat {a}}^{\dagger },\hat {a}$. Then, we multiply the kinetic equation from the right on vector  |I〉 and use (32)–(34) to convert the kinetic equation to the Schrödinger-like equation for the vector $\hspace{0.167em} \vert \rho (t)\rangle =\hat {\rho }(t)\hspace{0.167em} \vert I\rangle$:

$$\begin{eqnarray} &&\mathrm{i}\frac{\mathrm{d}}{\mathrm{d}t}\hspace{0.167em} \vert \rho (t)\rangle =L({\hat {a}}^{\dagger },\hat {a},{\tilde {a}}^{\dagger },\tilde {a})\hspace{0.167em} \vert \rho (t)\rangle , \end{eqnarray} \tag{ 35 }$$

where L is the Liouvillian which depends on both non-tilde and tilde superoperators. In particular, the Liouvillian for the Lindblad master equation (8) becomes

$$\begin{eqnarray} &&L=\hat {H}-\tilde {H}-\mathrm{i}\mathop{\sum }\limits _{b \alpha }{\Pi }_{b \alpha }, \end{eqnarray} \tag{ 36 }$$

where

$$\begin{eqnarray} {\Pi }_{b \alpha }&=&({\Gamma }_{b \alpha 1}-{\Gamma }_{b \alpha 2})({\hat {a}}_{b \alpha }^{\dagger }{\hat {a}}_{b \alpha }+{\tilde {a}}_{b \alpha }^{\dagger }{\tilde {a}}_{b \alpha })\nonumber\\ &&-~2\mathrm{i}({\Gamma }_{b \alpha 1}{\tilde {a}}_{b \alpha }{\hat {a}}_{b \alpha }+{\Gamma }_{b \alpha 2}{\tilde {a}}_{b \alpha }^{\dagger }{\hat {a}}_{b \alpha }^{\dagger })+2{\Gamma }_{b \alpha 2}. \end{eqnarray} \tag{ 37 }$$

In the derivation of (36) and (37) we took into account that $\hat {\rho }=\rho ({\hat {a}}^{\dagger },\hat {a})$ is a bosonic superoperator which commutes with all tilde superoperators. Due to the Lindblad dissipators, the Liouville superoperator (36) is non-Hermitian. In addition, as  |ρ〉 is tilde-invariant, the Liouvillian obeys the property $(L)^{\widetilde {}}=-L$.

Taking the time derivative of 〈I|ρ(t)〉 = 1 we find that 〈I| L = 0, i.e. 〈I|  the left zero-eigenvalue eigenstate of the Liouvillian superoperator. Since also 〈I|  is the vacuum for ${\hat {a}}_{j}^{\dagger }-\mathrm{i}{\tilde {a}}_{j}$ and ${\tilde {a}}_{j}^{\dagger }+\mathrm{i}{\hat {a}}_{j}$ superoperators, it is appropriate to call 〈I|  the left vacuum vector. For the electron transport problem we focus on nonequilibrium steady state where the current through the quantum dot is given by

$$\begin{eqnarray} &&\langle {J}_{\alpha }\rangle =\mathrm{Tr}({J}_{\alpha }{\rho }_{\infty })=\langle I\vert {\hat {J}}_{\alpha }\vert {\rho }_{\infty }\rangle . \end{eqnarray} \tag{ 38 }$$

Here ${\hat {J}}_{\alpha }$ is the current superoperator and the stationary, steady state solution of (35),  |ρ_∞〉, is the right zero-eigenvalue eigenstate (right vacuum vector) of the Liouville superoperator:

$$\begin{eqnarray} &&L \hspace{0.167em} \vert {\rho }_{\infty }\rangle =0. \end{eqnarray} \tag{ 39 }$$

In section 3, we show how one can find  |ρ_∞〉 perturbatively starting from the free-field approximation for the nonequilibrium density matrix.

Let us make the important remark on the notation used in the rest of the paper: only creation/annihilation operators written by letters a,d (such as, for example, a_bα and ${a}_{b \alpha }^{\dagger }$) are related to each other by the Hermitian conjugation; all other creation, c^†,b^†,γ^†, and annihilation c,b,γ, operators are 'canonically conjugated' to each other, i.e. for example, c^† does not mean (c)^† although cc^† ± c^†c = 1 (±—bosons/fermions). We will also use the same notation for the non-tilde superoperators ${\hat {a}}_{j}^{\dagger }$ and ${\hat {a}}_{j}$ as the ordinary operators ${a}_{j}^{\dagger }$ and a_j, bearing in mind that all operators acting in the Liouville–Fock space are superoperators.

We start by rewriting the Liouvillian (36) as

$$\begin{eqnarray} &&L={L}^{(0)}+{L}^{\prime}, \end{eqnarray} \tag{ 40 }$$

where L⁽⁰⁾ is the quadratic unperturbed part of L, and

$$\begin{eqnarray} &&{L}^{\prime}={H}_{\mathrm{S}}^{\prime}-{\tilde {H}}_{S}^{\prime} \end{eqnarray} \tag{ 41 }$$

is a perturbation. Then using the equation of motion method:

$$\begin{eqnarray} &&[{c}_{n}^{\dagger },{L}^{(0)}]=-{\Omega }_{n}{c}_{n}^{\dagger },\tqs [{c}_{n},{L}^{(0)}]={\Omega }_{n}{c}_{n}, \end{eqnarray} \tag{ 42 }$$

we exactly diagonalize[22] L⁽⁰⁾ in terms of nonequilibrium quasiparticle creation and annihilation operators:

$$\begin{eqnarray} &&{L}^{(0)}=\mathop{\sum }\limits _{n}({\Omega }_{n}c_{n}^{\dagger }{c}_{n}-{\Omega }_{n}^{\ast }{\tilde {c}}_{n}^{\dagger }{\tilde {c}}_{n}). \end{eqnarray} \tag{ 43 }$$

Here ${\tilde {c}}_{\sigma n}^{\dagger },{\tilde {c}}_{\sigma n}$ are obtained from ${c}_{\sigma n}^{\dagger },{c}_{\sigma n}$ by the tilde conjugation rules.

The nonequilibrium quasiparticle creation and annihilation operators are connected to ${a}^{\dagger },a,{\tilde {a}}^{\dagger },\tilde {a}$ by canonical (but not unitary) transformations:

$$\begin{eqnarray} \displaystyle {c}_{n}^{\dagger }&=&\mathop{\sum }\limits _{s}{\psi }_{n,s}b_{s}^{\dagger }+\mathop{\sum }\limits _{b \alpha }{\psi }_{n,b \alpha }b_{b \alpha }^{\dagger }, \\ \displaystyle {c}_{n}&=&\mathop{\sum }\limits _{s}({\psi }_{n,s}{b}_{s}+\mathrm{i}{\varphi }_{n,s}\tilde {b}_{s}^{\dagger })+\mathop{\sum }\limits _{b \alpha }({\psi }_{n,b \alpha }{b}_{b \alpha }+\mathrm{i}{\varphi }_{n,b \alpha }\tilde {b}_{b \alpha }^{\dagger }), \end{eqnarray} \tag{ 44 }$$

where

$$\begin{eqnarray*} &&{b}_{s}^{\dagger }={a}_{s}^{\dagger }-\mathrm{i}{\tilde {a}}_{s},\tqs {b}_{s}={a}_{s},\tqs {b}_{b \alpha }^{\dagger }={a}_{b \alpha }^{\dagger }-\mathrm{i}{\tilde {a}}_{b \alpha },\nonumber\\ &&{b}_{b \alpha }=(1-{f}_{b \alpha }){a}_{b \alpha }+\mathrm{i}{f}_{b \alpha }{\tilde {a}}_{b \alpha }^{\dagger }. \end{eqnarray*}$$

Nonequilibrium quasiparticle creation and annihilation operators obey the fermionic anti-commutation relations. In particular, from $\{ {c}_{n},{c}_{{n}^{\prime}}^{\dagger }\} ={\delta }_{n{n}^{\prime}}$ and $\{ {c}_{n},{\tilde {c}}_{{n}^{\prime}}\} =0$ we find the following orthonormality conditions for amplitudes:

$$\begin{eqnarray} &&\displaystyle \fl \mathop{\sum }\limits _{s}{\psi }_{n,\hspace{0.167em} s}{\psi }_{{n}^{\prime},\hspace{0.167em} s}+\mathop{\sum }\limits _{b \alpha }{\psi }_{n,\hspace{0.167em} b \alpha }{\psi }_{{n}^{\prime},\hspace{0.167em} b \alpha }={\delta }_{n{n}^{\prime}},\\ &&\displaystyle \fl \mathop{\sum }\limits _{s}({\psi }_{n,\hspace{0.167em} s}\varphi _{{n}^{\prime},\hspace{0.167em} s}^{\ast }-{\varphi }_{n,\hspace{0.167em} s}{\psi }_{{n}^{\prime},\hspace{0.167em} s}^{\ast })\\ && +~\mathop{\sum }\limits _{b \alpha }({\psi }_{n,\hspace{0.167em} b \alpha }\varphi _{{n}^{\prime},\hspace{0.167em} b \alpha }^{\ast }-{\varphi }_{n,\hspace{0.167em} b \alpha }{\psi }_{{n}^{\prime},\hspace{0.167em} b \alpha }^{\ast })=0. \end{eqnarray} \tag{ 45 }$$

By construction, 〈I| is the left vacuum for ${c}_{n}^{\dagger },{\tilde {c}}_{n}^{\dagger }$ operators. The vacuum for ${c}_{n},{\tilde {c}}_{n}$ operators, $\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle$, is automatically the zero-eigenvalue eigenstate of the unperturbed Liouvillian L⁽⁰⁾, i.e. it is the steady state density matrix in the zeroth-order approximation:

$$\begin{eqnarray} &&{L}^{(0)}\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle =0,\tqs \langle I\vert {\rho }_{\infty }^{(0)}\rangle =1. \end{eqnarray} \tag{ 46 }$$

In other words, the zeroth-order density matrix is the density matrix which does not contain nonequilibrium quasiparticle excitations.

Now we introduce the continuous real parameter λ, which will be set to unity at the end of the calculations:

$$\begin{eqnarray} &&L={L}^{(0)}+\lambda {L}^{\prime} \end{eqnarray} \tag{ 47 }$$

and expand the exact steady state density matrix in powers of λ:

$$\begin{eqnarray} &&\hspace{0.167em} \vert {\rho }_{\infty }\rangle =\mathop{\sum }\limits _{p=0}{\lambda }^{p}\hspace{0.167em} \vert \rho _{\infty }^{(p)}\rangle . \end{eqnarray} \tag{ 48 }$$

Substituting (48) into equation (39), we obtain the equation for the pth-order correction to the zeroth-order density matrix:

$$\begin{eqnarray} &&{L}_{0}\hspace{0.167em} \vert {\rho }_{\infty }^{(p)}\rangle =-{L}^{\prime}\hspace{0.167em} \vert {\rho }_{\infty }^{(p-1)}\rangle \end{eqnarray} \tag{ 49 }$$

or $\hspace{0.167em} \vert {\rho }_{\infty }^{(p)}\rangle =(-{L}_{0}^{-1}{L}^{\prime})^{p}\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle$. Here, L' is expressed in terms of the nonequilibrium quasiparticles. Thus, starting from  |ρ⁽⁰⁾〉 we can find any pth-order corrections to it. In addition, $\langle I\vert {\rho }_{\infty }^{(p)}\rangle =0$ for p ≥ 1 since $\hspace{0.167em} \vert {\rho }_{\infty }^{(p)}\rangle$ contains excited nonequilibrium quasiparticles.

To calculate the current through the quantum dot we express the current superoperator:

$$\begin{eqnarray} &&{J}_{\alpha }=-\mathrm{i}\mathop{\sum }\limits _{b s}{t}_{s b \alpha }(a_{b \alpha }^{\dagger }{a}_{s}-{a}_{s}^{\dagger }{a}_{b \alpha }) \end{eqnarray} \tag{ 50 }$$

in terms of nonequilibrium quasiparticle creation and annihilation operators and compute its expectation value with respect to 〈I|  and  |ρ_∞〉. As a result we get the following expansion:

$$\begin{eqnarray} &&{J}_{\alpha }=\mathop{\sum }\limits _{p=0}{\lambda }^{p}J_{\alpha }^{(p)}. \end{eqnarray} \tag{ 51 }$$

Here, ${J}_{\alpha }^{(0)}$ is the zeroth-order current for the system without interaction:

$$\begin{eqnarray} &&{J}_{\alpha }^{(0)}=-2\mathrm{Im}\mathop{\sum }\limits _{b s n}{t}_{s b \alpha }{\psi }_{n,b \alpha }{\varphi }_{n,s} \end{eqnarray} \tag{ 52 }$$

and ${J}_{\alpha }^{(p)}$ is the pth-order correction to it:

$$\begin{eqnarray} &&{J}_{\alpha }^{(p)}=-2\mathrm{Im}\mathop{\sum }\limits _{b s m n}{t}_{s b \alpha }\psi _{n,b \alpha }^{\ast }{\psi }_{m,s}{F}_{m n}^{(p)}, \end{eqnarray} \tag{ 53 }$$

where ${F}_{m n}^{(p)}$ is the expansion coefficient in

$$\begin{eqnarray} &&\hspace{0.167em} \vert {\rho }_{\infty }^{(p)}\rangle =\mathrm{i}\mathop{\sum }\limits _{m n}F_{m n}^{(p)}{c}_{m}^{\dagger }{\tilde {c}}_{n}^{\dagger }\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle +\cdots \end{eqnarray} \tag{ 54 }$$

and ${F}_{m n}^{(p)}=({F}_{n m}^{(p)})^{\ast }$ as follows from $\hspace{0.167em} \vert {\rho }_{\infty }^{(p)}\rangle ={\hspace{0.167em} \vert {\rho }_{\infty }^{(p)}\rangle }^{\widetilde {}}$. Thus, the problem of computing the pth-order correction to the unperturbed current is reduced to finding ${F}_{m n}^{(p)}$ by solving equation (49).

Using the same method we can obtain a perturbative expansion for the population of a quantum dot single-particle level:

$$\begin{eqnarray} &&{n}_{s}=\langle I\vert \hspace{0.167em} {a}_{s}^{\dagger }{a}_{s}\hspace{0.167em} \vert {\rho }_{\infty }\rangle =\mathop{\sum }\limits _{p=0}n_{s}^{(p)}, \end{eqnarray} \tag{ 55 }$$

where

$$\begin{eqnarray} &&\fl {n}_{s}^{(0)}=\mathop{\sum }\limits _{n}{\psi }_{n,s}{\varphi }_{n,s},\tqs n_{s}^{(p)}=-\mathop{\sum }\limits _{m n}{\psi }_{m,s}\psi _{n,s}^{\ast }{F}_{m n}^{(p)}.\nonumber\\ \end{eqnarray} \tag{ 56 }$$

The anti-commutation condition $\{ {b}_{s},{\tilde {b}}_{s}\} =0$ imposes the constraint on the amplitudes from which it follows that ${n}_{s}^{(0)}$ is a real number.

As the first application of the method we consider the Hamiltonian H_S which describes one electronic single-particle level coupled linearly to a vibration mode (phonon) of frequency ω₀ (the so-called local Holstein model):

$$\begin{eqnarray} &&{H}_{\mathrm{S}}={\varepsilon }_{0}{a}^{\dagger }a+{\omega }_{0}{d}^{\dagger }d+\kappa {a}^{\dagger }a({d}^{\dagger }+d). \end{eqnarray} \tag{ 57 }$$

For simplicity we assume that the tunneling matrix element in equation (5) is a real number t independent of indices α and b. The electron spin does not play any role here, so we suppress the spin index in the equations in this section. Replacing κ by λκ, we arrive at a perturbation expansion of the steady state density matrix  |ρ_∞〉 with respect to electron–vibronic coupling:

$$\begin{eqnarray} &&\hspace{0.167em} \vert {\rho }_{\infty }\rangle =\mathop{\sum }\limits _{p=0}{\lambda }^{p}\hspace{0.167em} \vert \rho _{\infty }^{(p)}\rangle . \end{eqnarray} \tag{ 58 }$$

To find the zeroth-order density matrix $\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle$, we diagonalize the fermionic part of L⁽⁰⁾. The resulting creation and annihilation operators have the form (44), and amplitudes ψ,φ satisfy the following system of equations:

$$\begin{eqnarray} && \displaystyle \hspace{0.167em} \hspace{0.167em} \hspace{0.167em} {\varepsilon }_{0}{\psi }_{n}-t\mathop{\sum }\limits _{b \alpha }{\psi }_{n,b \alpha }={\Omega }_{n}{\psi }_{n}\\ && \displaystyle \tqs {E}_{b \alpha }{\psi }_{n,b \alpha }-t{\psi }_{n}={\Omega }_{n}{\psi }_{n,b \alpha }, \end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray} && \displaystyle \hspace{0.167em} \hspace{0.167em} \hspace{0.167em} ({\varepsilon }_{0}-{\Omega }_{n}){\varphi }_{n}-t\mathop{\sum }\limits _{b \alpha }{\varphi }_{n,b \alpha }=t\mathop{\sum }\limits _{b \alpha }{f}_{b \alpha }{\psi }_{n,b \alpha }\\ && \displaystyle ({E}_{b \alpha }^{\ast }-{\Omega }_{n}){\varphi }_{n,b \alpha }-t{\varphi }_{n}=-t{f}_{b \alpha }{\psi }_{n}, \end{eqnarray} \tag{ 60 }$$

where E_bα = ε_bα − iγ_bα. The solution of eigenvalue problem (59) yields the spectrum of nonequilibrium quasiparticles, Ω_n and $-{\Omega }_{n}^{\ast }$, as well as ψ amplitudes which should be normalized according to equation (45). To find φ amplitudes we must solve linear equations (60).

Furthermore, let N_ω be the number of vibrational quanta with frequency ω₀ at temperature 1/β, i.e. N_ω = (exp(βω₀)−1)⁻¹. When κ = 0 the density matrix is factorized as $\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle =\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle _{f}\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle _{b}$:

$$\begin{eqnarray} &&\langle I\vert \hspace{0.167em} {d}^{\dagger }d \hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle ={N}_{\omega }. \end{eqnarray} \tag{ 61 }$$

It is convenient to introduce new phonon operators:

$$\begin{eqnarray} &&\gamma =(1+{N}_{\omega })d-{N}_{\omega }{\tilde {d}}^{\dagger },\tqs {\gamma }^{\dagger }={d}^{\dagger }-\tilde {d} \end{eqnarray} \tag{ 62 }$$

and their tilde conjugated partners, such that $\langle I\vert \hspace{0.167em} {\gamma }^{\dagger }=\langle I\vert \hspace{0.167em} {\tilde {\gamma }}^{\dagger }=0$ and $\gamma \hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle =\tilde {\gamma }\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle =0$. Then the quadratic part of the Liouvillian is diagonal in terms of introduced operators:

$$\begin{eqnarray} &&{L}^{(0)}=\mathop{\sum }\limits _{n}({\Omega }_{n}c_{n}^{\dagger }{c}_{n}-{\Omega }_{n}^{\ast }{\tilde {c}}_{n}^{\dagger }{\tilde {c}}_{n})+{\omega }_{0}({\gamma }^{\dagger }\gamma -{\tilde {\gamma }}^{\dagger }\tilde {\gamma }) \end{eqnarray} \tag{ 63 }$$

and the vacuum for ${c}_{n},{\tilde {c}}_{n},\gamma$, and $\tilde {\gamma }$ operators is the zeroth-order approximation for the density matrix, $\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle$. For the unperturbed current we have

$$\begin{eqnarray} &&{J}_{\alpha }^{(0)}=-2 t\mathrm{Im}\mathop{\sum }\limits _{b n}{\psi }_{n,b \alpha }{\varphi }_{n}, \end{eqnarray} \tag{ 64 }$$

while the pth-order correction is

$$\begin{eqnarray} &&{J}_{\alpha }^{(p)}=-2 t\mathrm{Im}\mathop{\sum }\limits _{b m n}\psi _{n,b \alpha }^{\ast }{\psi }_{m}{F}_{m n}^{(p)}. \end{eqnarray} \tag{ 65 }$$

To find ${F}_{m n}^{(p)}$ we rewrite the perturbative part of the Liouvillian in terms of operators c_n,γ, etc:

$$\begin{eqnarray} {L}^{\prime}&=&\mathop{\sum }\limits _{m n}\{ [L_{m n}^{(1)}{\gamma }^{\dagger }+{L}_{m n}^{(2)}{\tilde {\gamma }}^{\dagger }+{L}_{m n}^{(3)}(\gamma +\tilde {\gamma })]{c}_{m}^{\dagger }{c}_{n}-\mathrm{t.c.}\} \nonumber\\ &&-~\mathrm{i}\mathop{\sum }\limits _{m n}[L_{m n}^{(4)}{\gamma }^{\dagger }-({L}_{n m}^{(4)})^{\ast }{\tilde {\gamma }}^{\dagger }+{L}_{m n}^{(5)}(\gamma +\tilde {\gamma })]{c}_{m}^{\dagger }{\tilde {c}}_{n}^{\dagger }\nonumber\\ &&-~\mathrm{i}\mathop{\sum }\limits _{m n}L_{m n}^{(6)}({\gamma }^{\dagger }-{\tilde {\gamma }}^{\dagger }){c}_{m}{\tilde {c}}_{n}+\kappa {n}^{(0)}({\gamma }^{\dagger }-{\tilde {\gamma }}^{\dagger }), \end{eqnarray} \tag{ 66 }$$

where the coefficients L⁽ⁱ⁾ are

$$\begin{eqnarray} {L}_{m n}^{(1)}&=&\kappa [({\psi }_{m}-{\varphi }_{m})+{N}_{\omega }{\psi }_{m}]{\psi }_{n},\nonumber\\ \displaystyle {L}_{m n}^{(2)} &=& \kappa [{\varphi }_{m}+{N}_{\omega }{\psi }_{m}]{\psi }_{n},\tqs {L}_{m n}^{(3)}=\kappa {\psi }_{m}{\psi }_{n},\\ \displaystyle {L}_{m n}^{(4)} &=&\kappa [({\psi }_{m}-{\varphi }_{m}){\varphi }_{n}^{\ast }+{N}_{\omega }({\psi }_{m}{\varphi }_{n}^{\ast }-{\varphi }_{m}{\psi }_{n}^{\ast })]\\ \displaystyle {L}_{m n}^{(5)} &=& \kappa [{\psi }_{m}{\varphi }_{n}^{\ast }-{\varphi }_{m}{\psi }_{n}^{\ast }],\tqs {L}_{m n}^{(6)}={\psi }_{m}{\psi }_{n}^{\ast }, \end{eqnarray} \tag{ 67 }$$

and

$$\begin{eqnarray} &&{n}^{(0)}=\langle I\vert \hspace{0.167em} {a}^{\dagger }a \hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle =\mathop{\sum }\limits _{n}{\psi }_{n}{\varphi }_{n} \end{eqnarray} \tag{ 68 }$$

is an unperturbed electron level population. The notation 't.c.' in equation (66) means the tilde conjugation (i.e., ${c}_{m}^{\dagger }\rightarrow {\tilde {c}}_{m}^{\dagger }$, ${\gamma }^{\dagger }\rightarrow \tilde {\gamma }$, ${L}_{m n}^{(1)}\rightarrow ({L}_{m n}^{(1)})^{\ast }$, etc). Then, substituting equations (63) and (66) into (49) we obtain the following general expression for ${F}_{m n}^{(p)}$:

$$\begin{eqnarray} {F}_{m n}^{(p)}&=&-\frac{1}{{\Omega }_{m}-{\Omega }_{n}^{\ast }}\biggl \{ \mathop{\sum }\limits _{k}L_{m k}^{(3)}[{Z}_{k n}^{(p-1)}+({Z}_{n k}^{(p-1)})^{\ast }]\nonumber\\ &&-~\mathop{\sum }\limits _{k}(L_{n k}^{(3)})^{\ast }[({Z}_{k m}^{(p-1)})^{\ast }+{Z}_{m k}^{(p-1)}]-2{L}_{m n}^{(5)}{W}^{(p-1)}\biggr \} ,\nonumber\\ \end{eqnarray} \tag{ 69 }$$

where ${Z}_{m n}^{(p)}$ and W^(p) are coefficients in the expansion:

$$\begin{eqnarray} \hspace{0.167em} \vert {\rho }_{\infty }^{(p)}\rangle &=&\biggl \{ {W}^{(p)}({\gamma }^{\dagger }+{\tilde {\gamma }}^{\dagger })\nonumber\\ &&+~\mathrm{i}\mathop{\sum }\limits _{m n}\bigl [Z_{m n}^{(p)}{\gamma }^{\dagger }+({Z}_{n m}^{(p)})^{\ast }{\tilde {\gamma }}^{\dagger }\bigr ]{c}_{m}^{\dagger }{\tilde {c}}_{n}^{\dagger }+\cdots \biggr \} \hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle .\nonumber\\ \end{eqnarray} \tag{ 70 }$$

Thus, to find the pth-order correction to the current we need first to compute ${Z}_{m n}^{(p-1)}$ and W^(p−1). This can be done using the same method as used to find ${F}_{m n}^{(p)}$. As a result, ${Z}_{m n}^{(p)}$ and W^(p) are nonvanishing only for odd p. Therefore, only even powers of p contribute to the current expansion as it should be for the considered model. It is interesting to note that the term ${W}^{(p)}({\gamma }^{\dagger }+{\tilde {\gamma }}^{\dagger })\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle$ is associated with the momentum transfer from the electronic current to the quantum dot vibrational mode (current-induced translational motion of the dot) whereas ${Z}_{m n}^{(p)}{\gamma }^{\dagger }{c}_{m}^{\dagger }{\tilde {c}}_{n}^{\dagger }\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle$ and $({Z}_{n m}^{(p)})^{\ast }{\tilde {\gamma }}^{\dagger }{c}_{m}^{\dagger }{\tilde {c}}_{n}^{\dagger }\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle$ correspond to the current-induced heating and cooling processes, respectively.

As an example we give here explicit expressions for the first-order perturbation theory W⁽¹⁾ and ${Z}_{m n}^{(1)}$:

$$\begin{eqnarray} &&{W}^{(1)}=-\frac{{n}^{(0)}}{{\omega }_{0}},\tqs {Z}_{m n}^{(1)}=\frac{{L}_{m n}^{(4)}}{{\Omega }_{m}-{\Omega }_{n}^{\ast }+{\omega }_{0}}. \end{eqnarray} \tag{ 71 }$$

Combining equations (71) and (69), we find ${F}_{m n}^{(2)}$. Then inserting ${F}_{m n}^{(2)}$ into (65) we derive the second-order perturbation theory correction to ${J}_{\alpha }^{(0)}$. This correction consists of two parts: the first part is proportional to n⁽⁰⁾, so it is the Hartree term, while the remaining part is the Fock term. In section 3.4, we also verify these definitions by comparing Hartree and Fock terms obtained within the presented approach and the exact ones given by the NEGF formalism.

As a next example we consider electron transport through one spin-degenerate level with local Coulomb interaction:

$$\begin{eqnarray} &&{H}_{\mathrm{S}}={\varepsilon }_{0}\mathop{\sum }\limits _{\sigma }{n}_{\sigma }+U{n}_{\uparrow }{n}_{\downarrow }. \end{eqnarray} \tag{ 72 }$$

Here ${n}_{\sigma }={a}_{\sigma }^{\dagger }{a}_{\sigma }$ is the number operator for electrons with spin σ in the quantum dot. In what follows, we again assume the tunneling matrix element is independent of α,b as well as spin σ, i.e.

$$\begin{eqnarray} &&{H}_{\mathrm{SB}}=-t\mathop{\sum }\limits _{\sigma b \alpha }(a_{\sigma b \alpha }^{\dagger }{a}_{\sigma }+\mathrm{h.c.}). \end{eqnarray} \tag{ 73 }$$

We also assume that energy levels in the leads are spin-degenerate.

Since the quadratic part of the corresponding Liouvillian describes electron transport through noninteracting spin-up and spin-down levels, it is diagonalized by the same method as in the previous example. As a result we obtain

$$\begin{eqnarray} &&{L}^{(0)}=\mathop{\sum }\limits _{\sigma n}({\Omega }_{n}c_{\sigma n}^{\dagger }{c}_{\sigma n}-{\Omega }_{n}^{\ast }{\tilde {c}}_{\sigma n}^{\dagger }{\tilde {c}}_{\sigma n}). \end{eqnarray} \tag{ 74 }$$

The vacuum of c_σn and ${\tilde {c}}_{\sigma n}$ operators, $\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle$, is the density matrix in the zeroth-order perturbation theory and

$$\begin{eqnarray} &&{J}_{\alpha }^{(0)}=-4 t\mathrm{Im}\mathop{\sum }\limits _{b n}{\psi }_{n,b \alpha }{\varphi }_{n} \end{eqnarray} \tag{ 75 }$$

is the corresponding current.

To find the pth-order correction to (75):

$$\begin{eqnarray} &&{J}_{\alpha }^{(p)}=-4 t\mathrm{Im}\mathop{\sum }\limits _{b m n}\psi _{n,b \alpha }^{\ast }{\psi }_{m}{F}_{m n}^{(p)}, \end{eqnarray} \tag{ 76 }$$

we rewrite ${L}^{\prime}=U({n}_{\uparrow }{n}_{\downarrow }-{\tilde {n}}_{\uparrow }{\tilde {n}}_{\downarrow })$ in terms of nonequilibrium quasiparticles:

$$\begin{eqnarray} {L}^{\prime}&=&\mathop{\sum }\limits _{\sigma k l}\{ K_{k l}^{(1)}({c}_{\sigma k}^{\dagger }{c}_{\sigma l}-\mathrm{t.c.})+\mathrm{i}{K}_{k l}^{(2)}{c}_{\sigma k}^{\dagger }{\tilde {c}}_{\sigma l}^{\dagger }\} \nonumber\\ &&+~\mathop{\sum }\limits _{k l m n}\{ (L_{k l m n}^{(1)}{c}_{{k}_{\uparrow }}^{\dagger }{c}_{{l}_{\downarrow }}^{\dagger }{c}_{{m}_{\downarrow }}{c}_{{n}_{\uparrow }}-\mathrm{t.c.})+{L}_{k l m n}^{(2)}{c}_{k \uparrow }^{\dagger }{c}_{{l}_{\downarrow }}^{\dagger }{\tilde {c}}_{{m}_{\uparrow }}^{\dagger }{\tilde {c}}_{{n}_{\downarrow }}^{\dagger }\nonumber\\ &&+~{L}_{k l m n}^{(3)}({c}_{{k}_{\uparrow }}^{\dagger }{\tilde {c}}_{{l}_{\downarrow }}^{\dagger }{\tilde {c}}_{{m}_{\downarrow }}{c}_{{n}_{\uparrow }}+{c}_{{k}_{\downarrow }}^{\dagger }{\tilde {c}}_{l \uparrow }^{\dagger }{\tilde {c}}_{{m}_{\uparrow }}{c}_{{n}_{\downarrow }}\nonumber\\ &&-~{c}_{{k}_{\uparrow }}^{\dagger }{\tilde {c}}_{{l}_{\uparrow }}^{\dagger }{\tilde {c}}_{{m}_{\downarrow }}{c}_{{n}_{\downarrow }}-{c}_{{k}_{\downarrow }}^{\dagger }{\tilde {c}}_{{l}_{\downarrow }}^{\dagger }{\tilde {c}}_{{m}_{\uparrow }}{c}_{{n}_{\uparrow }})\nonumber\\ &&+~\mathrm{i}[{L}_{k l m n}^{(4)}({c}_{{k}_{\uparrow }}^{\dagger }{c}_{{l}_{\downarrow }}^{\dagger }{\tilde {c}}_{{m}_{\downarrow }}^{\dagger }{c}_{{n}_{\uparrow }}+{c}_{{k}_{\downarrow }}^{\dagger }{c}_{{l}_{\uparrow }}^{\dagger }{\tilde {c}}_{{m}_{\uparrow }}^{\dagger }{c}_{{n}_{\downarrow }})+\mathrm{t.c.}]\nonumber\\ &&+~\mathrm{i}[{L}_{k l m n}^{(5)}({c}_{{k}_{\uparrow }}^{\dagger }{\tilde {c}}_{{l}_{\downarrow }}{c}_{{m}_{\downarrow }}{c}_{{n}_{\uparrow }}+{c}_{{k}_{\downarrow }}^{\dagger }{\tilde {c}}_{{l}_{\uparrow }}{c}_{{m}_{\uparrow }}{c}_{{n}_{\downarrow }})+\mathrm{t.c.}]\} . \end{eqnarray} \tag{ 77 }$$

Here ${K}_{k l}^{(1)}$ and ${K}_{k l}^{(2)}$ are given by

$$\begin{eqnarray} \displaystyle {K}_{k l}^{(1)}&=&U{n}_{\sigma }^{(0)}{\psi }_{k}{\psi }_{l},\tqs {K}_{k l}^{(2)}=-U{n}_{\sigma }^{(0)}({\psi }_{k}{\varphi }_{l}^{\ast }-{\varphi }_{k}{\psi }_{l}^{\ast }), \\ \displaystyle {n}_{\sigma }^{(0)}&=&\langle I\vert \hspace{0.167em} {a}_{\sigma }^{\dagger }{a}_{\sigma }\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle =\mathop{\sum }\limits _{n}{\psi }_{n}{\varphi }_{n}, \end{eqnarray} \tag{ 78 }$$

while the coefficients ${L}_{k l m n}^{(i)}$ are listed in [23].

Now, substituting equations (74) and (77) into equation (49) we find the following general expression for ${F}_{m n}^{(p)}$:

$$\begin{eqnarray} {F}_{m n}^{(p)}&=&-\frac{1}{{\Omega }_{m}-{\Omega }_{n}^{\ast }}\biggl \{ {K}_{m n}^{(2)}{\delta }_{p 1}\nonumber\\ &&+~\mathop{\sum }\limits _{i}[K_{m i}^{(1)}{F}_{i n}^{(p-1)}\!-\!({K}_{n i}^{(1)}{F}_{i m}^{(p-1)})^{\ast }]\!-\!\mathop{\sum }\limits _{i j}\!L_{m n i j}^{(3)}{F}_{j i}^{(p-1)}\nonumber\\ &&-~\mathop{\sum }\limits _{i j k}[L_{m i j k}^{(5)}{G}_{k j n i}^{(p-1)}-({L}_{n i j k}^{(5)}{G}_{k j m i}^{(p-1)})]\biggr \} , \end{eqnarray} \tag{ 79 }$$

where δ_p1 is the Kronecker delta and ${G}_{k l m n}^{(p)}$ is a coefficient in the expansion:

$$\begin{eqnarray} &&\hspace{0.167em} \vert {\rho }_{\infty }^{(p)}\rangle =\biggl \{ \mathop{\sum }\limits _{k l m n}G_{k l m n}^{(1)}{c}_{\uparrow k}^{\dagger }{c}_{\downarrow l}^{\dagger }{\tilde {c}}_{\uparrow m}^{\dagger }{\tilde {c}}_{\downarrow n}^{\dagger }+\cdots \biggr \} \hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle . \end{eqnarray} \tag{ 80 }$$

In turn, the equation like (79) can be derived for ${G}_{k l m n}^{(p)}$.

The exact first-order perturbation theory correction to $\hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle$ is

$$\begin{eqnarray} \hspace{0.167em} \vert {\rho }_{\infty }^{(1)}\rangle &=&\biggl \{ \mathrm{i}\mathop{\sum }\limits _{\sigma m n}F_{m n}^{(1)}{c}_{\sigma m}^{\dagger }{\tilde {c}}_{\sigma n}^{\dagger }\nonumber\\ &&+~\mathop{\sum }\limits _{k l m n}G_{k l m n}^{(1)}{c}_{\uparrow k}^{\dagger }{c}_{\downarrow l}^{\dagger }{\tilde {c}}_{\uparrow m}^{\dagger }{\tilde {c}}_{\downarrow n}^{\dagger }\biggr \} \hspace{0.167em} \vert {\rho }_{\infty }^{(0)}\rangle , \end{eqnarray} \tag{ 81 }$$

where

$$\begin{eqnarray} \displaystyle {F}_{m n}^{(1)}&=& -\frac{{K}_{m n}^{(2)}}{{\Omega }_{m}-{\Omega }_{n}^{\ast }},\\ \displaystyle {G}_{k l m n}^{(1)}&=& -\frac{{L}_{k l m n}^{(2)}}{{\Omega }_{k}+{\Omega }_{l}-{\Omega }_{m}^{\ast }-{\Omega }_{n}^{\ast }}. \end{eqnarray} \tag{ 82 }$$

Inserting ${F}_{m n}^{(1)}$ into (76) we get the first-order perturbation theory correction ${J}_{\alpha }^{(1)}$ to the current (75). This correction is proportional to n⁽⁰⁾ and below we will show that it corresponds to the first-order Hartree term obtained with NEGF formalism.

Here we note that in [23] we applied perturbation theory to the Anderson model starting from the nonequilibrium Hartree–Fock approximation, i.e. L' was normal ordered and did not contain quadratic terms. Therefore, in [23] the mixture of two quasiparticle configurations to $\hspace{0.167em} \vert {\rho }_{\infty }^{(1)}\rangle$ vanished and the first-order perturbation theory correction to the current was zero.

To find the second-order correction to ${J}_{\alpha }^{(0)}$ we insert (82) into (79). This yields

$$\begin{eqnarray} {F}_{m n}^{(2)}&=&-\frac{1}{{\Omega }_{m}-{\Omega }_{n}^{\ast }}\biggl \{ \mathop{\sum }\limits _{i}[K_{m i}^{(1)}{F}_{i n}^{(1)}-({K}_{n i}^{(1)}{F}_{i m}^{(1)})^{\ast }]\nonumber\\ &&-~\mathop{\sum }\limits _{i j}L_{m n i j}^{(3)}{F}_{j i}^{(1)}-\mathop{\sum }\limits _{i j k}[L_{m i j k}^{(5)}{G}_{k j n i}^{(1)}-({L}_{n i j k}^{(5)}{G}_{k j m i}^{(1)})^{\ast }]\biggr \} .\nonumber\\ \end{eqnarray} \tag{ 83 }$$

Now, with the help of the obtained expression for ${F}_{m n}^{(2)}$ and equation (76) we get the second-order perturbation theory correction to ${J}_{\alpha }^{(0)}$.

Let us now compare the results obtained with the present approach with those calculated with the help of the Keldysh NEGF. For the case when the coupling to the left lead is proportional to the coupling to the right lead, the electron current through the quantum dot can be computed directly from the retarded dot Green's function, G^r(ω), using the well-known Meir–Wingreen formula [28]. For the considered models, assuming that the left and right leads are identical, Γ_L,R(ω) = 0.5Γ(ω), this formula takes the form

$$\begin{eqnarray} &&J=\frac{s}{2 \pi }\int \mathrm{d}\omega \hspace{0.167em} [{f}_{\mathrm{L}}(\omega )-{f}_{\mathrm{R}}(\omega )]\Gamma (\omega )\mathrm{Im}{G}^{\mathrm{r}}(\omega ). \end{eqnarray} \tag{ 84 }$$

Here s is the spin degeneracy of the considered models: s = 1 for the model with electron–vibration coupling and s = 2 for the model with electron–electron interaction. We will use the wide-band approximation for the electrode, so the imaginary part of the electrode self-energy, which is responsible for level broadening, is energy-independent, Γ(ω) = Γ, while its real part vanishes.

The retarded Green's function is the solution of the Dyson equation:

$$\begin{eqnarray} &&{G}^{\mathrm{r}}(\omega )={G}_{0}^{\mathrm{r}}(\omega )+{G}_{0}^{\mathrm{r}}(\omega ){\Sigma }^{\mathrm{r}}(\omega ){G}^{\mathrm{r}}(\omega ), \end{eqnarray} \tag{ 85 }$$

where ${G}_{0}^{\mathrm{r}}(\omega )=(\omega -{\varepsilon }_{0}+\mathrm{i}\Gamma )^{-1}$ is the noninteracting retarded Green's function and Σ^r(ω) is the retarded self-energy evaluated in the presence of electron–electron or electron–vibration interactions. Expanding Σ^r(ω) with respect to electron–electron or electron–vibration coupling, ${\Sigma }^{\mathrm{r}}(\omega )={\sum }_{p=1}{\lambda }^{p}{\Sigma }_{p}^{\mathrm{r}}(\omega )$, we obtain the perturbative expansion of G^r(ω) and consequently of the current:

$$\begin{eqnarray} J&=&\frac{s}{2 \pi }\int \mathrm{d}\omega \hspace{0.167em} [{f}_{\mathrm{L}}(\omega )-{f}_{\mathrm{R}}(\omega )]\Gamma (\omega )\nonumber\\ &&\times \ \mathrm{Im}\left [{G}_{0}^{\mathrm{r}}(\omega )+\mathop{\sum }\limits _{p=1}{\lambda }^{p}G_{p}^{\mathrm{r}}(\omega )\right ]=\mathop{\sum }\limits _{p=0}{\lambda }^{p}{J}^{(p)}. \end{eqnarray} \tag{ 86 }$$

Here J⁽⁰⁾ is the current through the system without interaction given by the standard Landauer formula.

In [22] we have shown that, for the current through a system without interaction, J⁽⁰⁾, the exact agreement between NEGF and the kinetic equation approach can be achieved by increasing the density of states in the buffer zones. Below we show that this is also true for correlated electronic systems.

In what follows, in the calculations based on the kinetic equation we will assume that N single-particle levels in each buffer zone are evenly distributed within the energy bandwidth [E_min;E_max] = [ − 10,10]. This bandwidth is much larger than any energy parameter in the system, so it corresponds to the wide-band approximation used in NEGF calculations. The tunneling coupling strength t is computed from Γ = 2πηt², where η = N/(E_max − E_min) is the density of states in the buffer zone. We note here that the main approximation in the derivation of the Lindblad master equation (8) is that the single-particle states in the buffer zone propagate in time as free states (7). It is evident from equation (7) that the larger the buffer zone, i.e. the larger the density of states η, the better this approximation. This will also be explicitly demonstrated in the numerical calculations below. The parameter γ in the Lindblad operators is chosen to be γ = 2Δε, where ε is the energy spacing between the energy levels in the buffer zone. In addition, although it is not necessary, we use a symmetric applied voltage, μ_L,R =± 0.5 V , and the temperature of the electrodes is zero, T = 0.

Initially, we consider the system with electron–vibronic interaction and compare the second-order correction to the current obtained in section 3.2 with that calculated using NEGF formalism (86). We use the following model parameters of the Hamiltonian (57): κ = 1.0,ω₀ = 1.0.

In NEGF formalism the second-order correction to the current arises from the retarded self-energy ${\Sigma }_{2}^{\mathrm{r}}$ which contains contributions from Hartree and Fock diagrams, ${\Sigma }_{2}^{\mathrm{r}}={\Sigma }_{\mathrm{H}}^{\mathrm{r}}+{\Sigma }_{\mathrm{F}}^{\mathrm{r}}$. The Hartree self-energy is [8]

$$\begin{eqnarray} &&{\Sigma }_{\mathrm{H}}^{\mathrm{r}}(\omega )=-\frac{2{\kappa }^{2}}{{\omega }_{0}}{n}^{(0)}, \end{eqnarray} \tag{ 87 }$$

where n⁽⁰⁾ is the electron level population in the zero-order approximation:

$$\begin{eqnarray} &&{n}^{(0)}=\frac{\Gamma }{2 \pi }\int \mathrm{d}\omega \frac{{f}_{\mathrm{L}}(\omega )+{f}_{\mathrm{R}}(\omega )}{(\omega -{\varepsilon }_{0})^{2}+{\Gamma }^{2}}. \end{eqnarray} \tag{ 88 }$$

The expression for the Fock self-energy is more complicated and can be found elsewhere (see, for example, [29]).

In figures 2 and 3 we compare Hartree and Fock second-order corrections to the current obtained within our approach with different size N of buffer zone and the exact ones. The corrections are shown as functions of the level energy, ε₀, for two values of the applied voltage V and broadening Γ. It is evident from the figures that the difference between exact and Lindblad-equation-based results become negligible as we increase the lead density of states in the buffer zone. The reason is that increasing the number of single-particle states in the buffer zones we make approximation (iii), under which Lindblad master equation (8) was derived, more justified. The deviation of the results obtained from the Lindblad kinetic equation and NEGF becomes smaller at the larger applied voltage or Γ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Now we compare first-order corrections to the current for the Anderson model. We put U = 1.0 for the strength of the Coulomb interaction. Within the NEGF formalism the first-order correction is solely due to the Hartree diagram and it is

$$\begin{eqnarray} &&{J}^{(1)}=4{\Gamma }^{2}U{n}^{(0)}\int \frac{\mathrm{d}\omega }{2 \pi }\frac{({f}_{\mathrm{L}}(\omega )-{f}_{\mathrm{R}}(\omega ))(\omega -{\varepsilon }_{0})}{((\omega -{\varepsilon }_{0})^{2}+{\Gamma }^{2})^{2}}, \end{eqnarray} \tag{ 89 }$$

where the population n⁽⁰⁾ is given by equation (88).

The results of numerical calculations are shown in figure 4 for different values of Γ and applied voltage V. As we can see the results of the Lindblad equation approach and converge to the exact results with increasing value of N and the convergence is faster for larger values of applied voltage and Γ.

Zoom In Zoom Out Reset image size

In figure 5 we show the current through the Anderson model computed by means of the Lindblad equation by taking into account the first- and second-order corrections. We take N = 1600, so the obtained results correspond to NEGF ones. As we can see from the figure, the first- and second-order contributions shift the maximum of the current towards the symmetric point ε₀ =− 0.5U. The first-order correction increase the maximum current, while the second-order correction acts in the opposite direction. We also see from figure 5 that for a given U the relative value of the first- and second-order corrections show little dependence on the applied voltage V. In contrast, in [23] we have observed that nonequilibrium post-Hartree–Fock electronic correlations play the important role at larger applied voltages and, as a result, the second-order correction to the current become more pronounced with increasing V. This is due to the difference in the structure and spectrum of nonequilibrium quasiparticles. The quasiparticle spectrum, both ψ and φ amplitudes, depend on the voltage in the post-Hartree–Fock perturbation theory [23], whereas in the present work the voltage enters only into φ amplitudes of the nonequilibrium quasiparticles through Fermi–Dirac occupation numbers of the buffer states.

Zoom In Zoom Out Reset image size