# Density-operator approaches to transport through interacting quantum dots: simplifications in fourth order perturbation theory

###### Abstract

Various theoretical methods address transport effects in quantum dots beyond single-electron tunneling, while accounting for the strong interactions in such systems. In this paper we report a detailed comparison between three prominent approaches to quantum transport: the fourth order Bloch-Redfield quantum master equation (BR), the real-time diagrammatic technique (RT) and the scattering rate approach based on the T-matrix (TM). Central to the BR and RT is the generalized master equation for the reduced density matrix. We demonstrate the exact equivalence of these two techniques. By accounting for coherences (non-diagonal elements of the density matrix) between non-secular states, we show how contributions to the transport kernels can be grouped in a physically meaningful way. This not only significantly reduces the numerical cost of evaluating the kernels, but also yields expressions similar to those obtained in the TM approach, allowing for a detailed comparison. However, in the TM approach an ad-hoc regularization procedure is required to cure spurious divergences in the expressions for the transition rates in the stationary (zero-frequency) limit. We show that these problems derive from incomplete cancellation of reducible contributions and do not occur in the BR and RT techniques, resulting in well-behaved expressions in the latter two cases. Additionally, we show that a standard regularization procedure of the TM rates employed in the literature does not correctly reproduce the BR and RT expressions. All the results apply to general quantum dot models and we present explicit rules for the simplified calculation of the zero-frequency kernels. Although we focus on fourth order perturbation theory only, the results and implications generalize to higher orders. We illustrate our findings for the single impurity Anderson model with finite Coulomb interaction in a magnetic field.

###### pacs:

73.23.Hk, 73.63.-b, 73.63.Kv, 73.40.Gk## I Introduction

The experimental progress in fabrication of ultrasmall electrical devices Tarucha et al. (1996); Tans et al. (1997); Martel et al. (1998); Tans et al. (1998); Collier et al. (1999); Bachtold et al. (2001); Huang et al. (2001) has made quantum dots one of the standard components in fundamental research and application oriented nanostructures. Whereas high-resolution transport measurements in the low temperature regime have reached a high degree of sophistication and reveal data dominated by complex many-body phenomena Franceschi et al. (2001); Liang et al. (2002); Park et al. (2002); van der Wiel et al. (2003); Yu and Natelson (2004); Schleser et al. (2005); Jarillo-Herrero et al. (2005); Sapmaz et al. (2006); Osorio et al. (2007); Parks et al. (2007); Hauptmann et al. (2008), theoretical methods are still struggling to describe these, mainly due to competing influences of strong local interactions and quantum fluctuations Schrieffer and Wolff (1966); Anderson (1970); Anderson et al. (1970); König et al. (1997, 1998); Paaske et al. (2004); Pustilnik and Glazman (2004); Thorwart et al. (2005); Schoeller and Reininghaus (2009); Koch et al. (2006); Holm et al. (2008); Leijnse and Wegewijs (2008); Leijnse et al. (2009); Begemann et al. (2010); Andergassen et al. (2010).

A common setup for transport studies is drawn in Fig. 1(a): the number of electrons on the device is controlled capacitively via a gate voltage , a difference in the electro-chemical potentials of the leads is created by a bias voltage .
The measured quantity is the current or the differential conductance of the whole circuit, which is usually represented in stability diagrams, where the changes in current vs. and are color-coded [Fig. 1(b)].
The observed phenomena strongly depend on the strength of the coupling of the nanodevice to the electronic reservoirs.
In the limit of extremely weak coupling, current at low bias is completely blocked in wide ranges of the gate voltage, showing up as so-called Coulomb diamonds [e.g. Fig. 1(b), central region].
Outside these regions of Coulomb blockade only single electrons can be transferred sequentially onto or out of
the dot Ingold and Nazarov (1992); Kouwenhoven et al. (1997), a process called single-electron tunneling (SET).
For simple systems (e.g. without orbital degeneracies) in this regime rate equations H.J. Kreuzer (1981)
are the standard technique to calculate the occupations of the dot states, the current and other transport quantities Bruus and Flensberg (2004).
The transition rates are calculated by Fermi’s Golden Rule, i.e., leading order perturbation theory in the tunneling.
For more complex quantum dots with degenerate orbitals Wunsch et al. (2005); Harbola et al. (2006); Mayrhofer and Grifoni (2006); Donarini et al. (2006); Begemann et al. (2008); Donarini et al. (2009); Schultz and von Oppen (2009); Schultz () and/or non-collinear magnetic electrodesBraun et al. (2004); Weymann and Barnas (2007); Koller et al. (2007); Hornberger et al. (2008), coherences,
i.e., non-diagonal density matrix elements, give crucial contributions to the
transport quantities and cannot be neglected.
These are typical situations in molecular electronics and spintronics Reckermann et al. (2009).

Since the transparency of the contacts is a matter of the material choice as well as of fortune, on the way to low ohmic contacts, intermediate coupling strengths are often observed, allowing for coherent tunneling of multiple electrons D. V. Averin and Y. V. Nazarov (1992).
Also, there is the possibility to design structures with tunable tunnel barriers, such that different coupling regimes can be systematically accessed, allowing for more detailed spectroscopic information to be extracted Franceschi et al. (2001).
Therefore, electron transport theory must go beyond lowest order perturbation theory in the tunneling,
while including many charge states, their complex excitations and their quantum coherence.
In recent years, several advanced approaches that address higher-order effects have been developed based on iterative real-time path-integral methods S. Weiss et al. (2008), scattering-states Mehta and Andrei (2006) combined with quantum-monte Carlo Han and Heary (2007) or numerical R. Bulla et al. (2008); Anders (2008) or analytical renormalization group methods Paaske et al. (2004); Schoeller and Reininghaus (2009).
Although these new methods are promising, the standard generalized master equation (GME) approach
still offers several advantages.
The GME describes the reduced density matrix of the quantum dot with transport kernels that are calculated perturbatively.
The GME can be derived using various methods Timm (2008): the Nakajima-Zwanzig Nakajima (1958); Zwanzig (1960)
projection operator technique Breuer and Petruccione (2002); Fick and Sauermann (1990); Kühne and Reinecker (1978), the real-time diagrammatic technique Schoeller and Schön (1994); König et al. (1997); Schoeller (1997) (RT) and the Bloch-Redfield approach Wangsness and Bloch (1953); Bloch (1957); Redfield (1965) (BR).
Evaluating the kernels up to fourth order in the tunneling Hamiltonian (next-to-leading order),
one can account for all processes involving coherent tunneling of one or two electrons.
These corrections to SET can be calculated either analytically for simple models
or, in complex cases, in a numerically efficient way.
In this case the GME is clearly limited to moderate values of the tunnel coupling as compared to temperature.
However, it has the benefit of non-perturbatively treating both the interactions on the dot as well as the non-equilibrium conditions imposed by the
bias voltage.
It can therefore provide crucial physical insights into measurements of non-linear transport through complex quantum dots,
see e.g. Ref. Hüttel et al. (2009); Zyazin et al. (2010).
More generally, even higher corrections can be explicitly formulated in any order of the tunneling by systematic diagram rules.
Since the explicit form of the kernel is known in this way,
a renormalization-group theory for its calculation can be formulated as well,
allowing the non-equilibrium low-temperature regime to be addressed Schoeller (2009), including the Kondo effect Schoeller and Reininghaus (2009).
A class of contributions beyond fourth order can also be included by expanding an equation of motion for the density
matrix Nyvold Pedersen and Wacker (2005, 2010).

This paper focuses on the BR and RT formulations of the GME approach applied to transport in the fourth order of perturbation theory and addresses crucial technical matters and simplifications relevant for the description of complex quantum dots. Several important concrete issues have motivated this work:

(i) There is an ongoing discussion about the validity and equivalence of approaches which can become obscured by the complexity of the expressions involved when discussing complex quantum dots. Clearly, the general form of the quantum master equation is well known since several decades. Still, a much debated issue is the actual task of systematically calculating higher order corrections to the transport kernels occurring in this equation for general complex quantum dot models. This paper emphasizes that the BR and RT techniques are one-to-one equivalent. In contrast, the scattering rate approach based on the generalized Fermi’s golden rule and T-matrix (TM), as formulated in the literature, differs from these two techniques Timm (2008). Although it also relies on a fourth order perturbative calculation, the results do not coincide in general for identical models. The reason is that the objects calculated perturbatively, i.e. the T-matrix and the time evolution kernel, respectively, are different objects whose relation needs to be clarified. In particular, the divergences occurring in the TM method are intrinsic to the method and not to the problem. We show that these go back to a lack of cancellation of divergent, reducible contributions to the transport kernels and that the regularizations proposed in the literature cannot reproduce the exact fourth order kernel. We quantitatively demonstrate the resulting deviations from the correct GME (BR or RT) result for the example of a single impurity Anderson model in magnetic field and analytically show how the divergent TM expressions are automatically regularized in the GME approaches. The GME approaches consistently account for all contributions to the perturbation expansion of the transport kernels in a given order. The importance of this was recently highlighted for the well studied non-equilibrium Anderson model, which was found to exhibit a previously unnoticed resonance due to coherent tunneling of electron pairs Leijnse et al. (2009).

(ii) The importance of non-diagonal elements in lowest order calculations involving degenerate states has long been recognized (“secular contributions”), and continues to attract attention in the context of transport. Only recently, the importance of non-secular terms (coherences between non-degenerate states) was found to be crucial Leijnse and Wegewijs (2008) for fourth order tunnel effects. We generalize the discussion in Ref. Leijnse and Wegewijs (2008) and show how these non-secular corrections can efficiently be included into effective fourth order transport kernels through certain reducible diagrams.

(iii) Explicit expressions for the fourth order transport kernels for a very general class of quantum dots were derived in Ref. Leijnse and Wegewijs (2008). However, the numerical cost of evaluating these expressions limits their applicability to systems where a relatively small number of many-body excitations () has to be accounted for. Here we show how contributions to the effective kernels can be grouped, making generally valid cancellations explicit and resulting in fewer and simpler terms in the perturbation expansion. From direct comparison between numerical implementations of the expressions in Ref. Leijnse and Wegewijs (2008) and of our new ”grouped” expressions, we find the latter to be between 10 and 20 times faster, without introducing any additional approximation. This allows the treatment of more complicated and realistic quantum dot models. The direct gain due to the reduction of the number of expressions amounts to a decrease of the computation time by a factor of 4. However, the grouping structure can be exploited further to make the numerical implementation more efficient, leading to the additional speed-up. Moreover, the grouping gives a basis upon which an explicit connection to TM expressions can be revealed.

As will become clear in the course of this paper, our newly found grouping intimately connects and enlightens the above three issues, which warrants our systematic and extended discussion. The key ideas presented can be applied to analyze higher order contributions as well as to similar perturbation and renormalization group calculations for other classes of problems.

The structure of the paper is as follows. In Sec. II we discuss the model Hamiltonian of the setup Fig. 1(a) and some pertinent notation. We then introduce the reduced density matrix (RDM) describing the quantum dot as part of the whole system and the generalized master (or kinetic) equation (GME) which describes its time-evolution. We summarize its general properties and the common ground of the discussed approaches. We then turn to the derivation of the generalized master equation using the BR and RT technique. The crucial role played by time-ordering, irreducibility of contributions and analytic properties (lack of spurious divergences) is emphasized. The derivations are given as compactly as possible, because there exists a broad formal study on different master-equation approaches by Timm Timm (2008), who also showed their equivalence. In contrast to his work, our comparison continues in Sec. III on a more explicit level with a mapping between terms arising from the BR and RT diagrams. In Sec. IV we demonstrate how the non-secular contributions of coherences lead to important corrections in the fourth order transport rates. Based on this, Sec. V introduces a grouping of contributions to the transport kernels, yielding significant simplifications due to partial cancellations, which is followed by an analysis of how the groups of diagrams contribute to fourth order physical transport processes (cotunneling, pair tunneling, and level renormalization and broadening). In Sec. VI the derivation of the TM is reformulated. The long-standing problem of a precise comparison with the RT technique in the context of transport theory and the origin of divergences in the TM approach is solved for the general case. The theoretical discussion in the last two sections, V and VI, is illustrated by the tangible application to a single impurity Anderson model in magnetic field.

## Ii Model and generalized master equation

The standard model for a quantum dot system coupled to contacts reads

(1) |

The Hamiltonian

(2) |

models the reservoirs, i.e., the source and the drain contact.
The operator () creates (annihilates) an electron in a
state with energy in the source () or drain () contact,
where denotes the spin projection.
The bias voltage shifts the electro-chemical potentials of the source and drain leads such that , where is the electron charge.

The coupling between the quantum dot and the leads is described by the tunnel Hamiltonian

(3) |

where () creates (annihilates) an electron in the single particle state on the dot. The single-particle amplitude for tunneling from an orbital state in lead to an orbital state on the dot is assumed to be independent of spin. Finally, the dot is described by the Hamiltonian

(4) |

where is a many-body eigenstate of the dot with energy . The precise dependence of these energies on the applied voltages arising from capacitive effects (see e.g., Ref. Kouwenhoven et al. (1997); von Delft and Ralph (2001)) is irrelevant for the following discussion. Typically, the gate voltage dependence is linear, , where is the number of electrons for state .

The diagonalized many-body Hamiltonian (4) together with the tunnel matrix elements (TMEs) between all the many-body eigenstates ,

(5a) | |||||

(5b) |

form the crucial input to the GME transport theory, which thereby incorporates local interaction effects non-perturbatively. Here is only nonzero if the state differs from in its electron number by . Typically, the TMEs can be assumed independent of , but this is not a prerequisite for the results of this work.

We focus on the regime of weak tunnel coupling and thus split in a free part and a perturbation . The condition for weak coupling is that the broadening induced by tunneling processes from and to lead is small compared to the thermal energy, i.e., , where and are respectively the Planck and Boltzmann constants. The broadening is defined by

(6) |

where for the contexts in which we need it here (i.e., as a measure of order of magnitude), both the orbital () dependence and the frequency () dependence are neglected. Notice that the analytical expressions we derive in this work are a priori not subject to these restrictions. For simplicity, we will further denote contributions of the order of in general by . In this paper we will go beyond lowest order in , allowing the regime of intermediate coupling to be addressed.

Throughout the paper, we will make use of the Liouville (superoperator) notation, in addition to standard Hamiltonian operator expressions. The former is merely an efficient bookkeeping tool on a general level, whereas the latter may be more useful when evaluating matrix elements. For example, for the dot Hamiltonian the abbreviation

(7) |

defines the action of the Liouville superoperator in the Schrödinger picture on an arbitrary operator . It generates the time-evolution through (Baker-Campbell-Hausdorff formula). Analogous expressions hold for the other Hamiltonians, and in particular we will need

(8) |

where is an operator in the interaction picture. Notice that the time evolution of the superoperator in the interaction picture is thus

(9) |

### ii.1 Generalized master equation and steady state

The object of interest here is the reduced density matrix Blum (1996) (RDM)

(10) |

It describes the state of the quantum dot incorporating the presence of the leads, which are traced out of the total density matrix , as prescribed by . Once is known, the expectation value of any observable can be calculated, as discussed below. When the interaction is switched on at time , the total density matrix is the direct product of the (arbitrary) initial state of the quantum dot and the equilibrium state of the leads,

(11) |

with . After the interaction is switched on, i.e. for times , correlations, which are of the order of the tunnel coupling Blum (1996), build up between leads and quantum dot, causing to deviate from the factorized form:

(12) |

We emphasize that it is crucial to include in a kinetic equation for , the correlations between leads and quantum dot consistently beyond linear order in , if one is interested in going beyond lowest order. As we will see in Secs. II.2 and II.3, the RT approach incorporates them automatically by directly integrating out the leads for times , while for the BR one explicitly solves for the deviation from the factorized state. Both the BR and the RT technique lead to the generalized quantum master equation (or kinetic equation), describing the time evolution of the RDM

(13) |

Here, the first term accounts for the time evolution due to the local dynamics of the quantum dot. In the second time non-local term, the time evolution kernel is a superoperator acting on the density operator. Convoluted in time with , it gives that part of the time evolution which is generated by the tunneling. We note that this form of the GME is dictated by the linearity of the Liouville equation and the partial trace operation.

The kernel to fourth order formally reads

(14) | ||||

Below we will show how BR and RT (as well as the Nakajima-Zwanzig projection operator approach, see App. A), when consistently applied, lead to this result. The central topic of this paper is the explicit evaluation of Eq. (14) and the significant simplifications which can be achieved in the steady state limit . In this limit, Eq. (13) becomes

(15) |

where and denotes an imaginary infinitesimal, and

(16) |

is the Laplace transform of the time evolution kernel. Taking matrix elements with respect to the many-body eigenstates of the dot Hamiltonian, , we obtain from Eq. (15) a set of linear coupled equations for all states of the RDM:

(17) |

Here, the matrix elements of (or any other superoperator) are defined by

(18) |

where we use square brackets to make clear that the kernel superoperator must first act on , and then the matrix elements of the resulting operator are taken. Each diagonal element of the RDM equals the probability of finding the system in a certain state. Thus, the normalization condition

(19) |

must be fulfilled. The restriction (19) allows the system of linear equations obtained from Eq. (17) to be solved, since without it they are under-determined due to the sum-rule

(20) |

Physically, this guarantees that gain and loss of probability are balanced in the stationary state.

The expectation value of any non-local observable can be expressed in a form similar to Eq. (13). In particular, we can write the particle current flowing out of lead (i.e. the number of electrons leaving lead per unit time) as

(21) |

where is the kernel associated with the current operator

(22) |

with being the number operator in lead .

Taking the steady state limit of Eq. (21), the stationary current is given by the zero-frequency component
of the Laplace transform of the current kernel, traced over in product with
the stationary density matrix :

(23) |

The current kernel can be calculated by simple modification of the time-evolution kernel as discussed below in subsections II.2 and II.3 explicitly.
We will now address the derivation of Eq. (13) and of its kernel (14)
up to fourth order in the tunnel coupling.
We focus here on two approaches, an iterative procedure in the time
domain Wangsness and Bloch (1953); Bloch (1957); Redfield (1965), referred to as Bloch-Redfield approach
(BR), and the RT approach Schoeller and Schön (1994); König et al. (1997); Schoeller (1997)
(RT).
The projection operator technique of NakajimaNakajima (1958) and Zwanzig Zwanzig (1960),
which has been explained and used in many works Breuer and Petruccione (2002); Fick and Sauermann (1990); Kühne and Reinecker (1978); Jang et al. (2001),
is closely related and equivalent to the BR approach
and is discussed for completeness in App. A.

The derivation of the kinetic equation requires no other
ingredient than the Liouville equation for the total density matrix Blum (1996):

(24) |

For the purposes of this paper, it is most convenient to work in the time-domain and use the interaction picture. In addition, we make use of the property of the particular bi-linear coupling of Eq. (3) considered here, that the lead-average of an odd number of interactions vanishes due to the odd number of lead field operators in .

### ii.2 Bloch-Redfield approach

The Bloch-Redfield approach Wangsness and Bloch (1953); Bloch (1957); Redfield (1965) is usually favored to derive the second order quantum master equationBlum (1996). Basically, one integrates Eq. (24) and reinserts it back into its differential form to get

(25) |

We now extend this to fourth order Laird et al. (1990) by repeating the iteration steps: we transform Eq. (25) to an integral equation,

(26) |

which is once more reinserted into Eq. (24). After integration one arrives at

(27) |

We reinsert Eq. (27) back into the Liouville equation (24) and perform the trace over the leads in order to obtain the RDM. Thereby, terms involving in total an odd number of lead operators vanish. Due to the relations and with Eq. (12) we obtain:

(28) |

The second order contribution in Eq. (28) contains instead of and lacks the convoluted form which the fourth order term has: thus in the stationary limit the initial state does not seem to drop out. If one naively were to neglect this difference and set the fourth order kernel would contain spurious divergences (see Sec. VI). Instead, one has to account for the correlations between dot and reservoirs at times up to order as expressed by Eq. (26). Taking the trace over the leads this equation gives:

(29) |

This shows also that Eq. (28) through still contains higher-order terms. Consistently neglecting these in Eq. (28), we can eliminate the dependence on the initial condition from Eq. (28) and thus arrive at the GME in the interaction picture

(30) |

with the time-evolution kernel defined by

(31) |

Transforming the RDM to the Schrödinger picture by

(32) |

and with the Liouville operators according to Eq. (9),
we obtain the generalized master equation (13),
and arrive at the expression (14) for the kernel.

The current, Eq. (21), is given by

(33a) | ||||

(33b) |

where the current kernel in the interaction picture is given by Eq. (31) with replaced by . In deriving this, as for the density matrix, one must take care to keep the time-ordered structure: since the current operator is, as , linear in the lead operators [cf. Eq. (22)], we obtain Eq. (33b), correct up to fourth order, by inserting the third order iteration for (Eq. (27)) into Eq. (33a). Under the trace, the first and third contribution are zero since they contain an odd number of interactions. Next, as before, in the second contribution of Eq. (27), has to be eliminated using Eq. (29), thereby generating a fourth order correction term. Finally, in the fourth term one must consistently keep , i.e. only here one can neglect the deviation from the factorized form.

### ii.3 Real-time diagrammatic technique

The real-time approach has been discussed on a general level in many works Schoeller and Schön (1994); König et al. (1997); Schoeller (1997). Therefore the aim of this section is to recall how one efficiently arrives at the kinetic equation and its kernel by exploiting Wick’s theorem from the outset. We start from the Liouville equation (24) for the full system in the interaction picture and formally integrate it:

(34) |

where is the time-ordering superoperator. Using and defining the superoperator

(35) |

the time-evolution of the reduced density matrix can formally be written as:

(36) |

Expanding the time-ordered exponential superoperator, the trace in Eq. (35) can be explicitly evaluated term by term by Wick’s theorem: the trace over each string of reservoir field operators becomes a product of pair contractions, indicated in the following by contraction lines. For our purposes here, one can simply formally consider the Liouvillians to be contracted (meaning their reservoir part), see Eq. (38) below. We can then decompose into a reducible and an irreducible part, depending on whether or not the contractions separate into disconnected blocks. Collecting all irreducible parts into the kernel one obtains in the standard way a Dyson equation:

(37) |

It relates the full propagator to the free propagator, which equals unity in the interaction picture, and to the irreducible kernel .
Applying the Dyson equation to and taking the time derivative, one arrives at the kinetic equation in the interaction picture, Eq. (30). Transformed back to the Schrödinger picture we obtain the kinetic equation (13).

We have now obtained the kernel formally as the sum of irreducible contributions to the time-evolution superoperators of different orders in the tunneling, which are written down to fourth order:

(38) |

where .

The expectation value of the current (or of any operator) is obtained in a similar fashion. We first introduce a superoperator , which is an anti-commutator in contrast to the superoperators governing the time-evolution. The current is then expressed as

(39) |

where we introduced a current propagator

(40) |

which differs from the propagator for the reduced density operator only by the current vertex at the final time. Collecting all parts of which are irreducibly connected to the latter vertex, one readily verifies that the remaining irreducible parts at earlier times are those contained in the propagator :

(41) |

The current kernel is given formally as the sum of irreducible contributions to the time-evolution superoperators of different orders in the tunneling with the leftmost interaction vertex replaced by the current superoperator. Applying this equation to and taking the trace over the dot, one arrives at the expression for the current in terms of the new kernel and the reduced density matrix in the interaction picture, Eq. (33b). Transformed back to the Schrödinger picture we obtain Eq. (21).

### ii.4 Comparison of the approaches

For the comparison between the BR and RT approaches, it is most useful to consult Eq. (31) and Eq. (38). The second order terms, contained in both equations in the first line, obviously match. The equivalence of the fourth-order terms is more indirect: in the BR approach, the first term of Eq. (31) is evaluated using Wicks’ theorem by building all possible contractions, including the reducible ones (contraction between vertices at times and as well as and ). The latter are then canceled by the second term. Precisely the same happens in the projection operator approach [cf. Eq. (77)]. The above conclusions hold in fact for any order of perturbation theory as shown in Timm (2008). In contrast, the RT approach avoids the inclusion and subsequent cancellation of reducible parts which rapidly grow in number with the perturbation order.

We emphasize that there is one unique correct fourth order (time non-local) generalized master equation, in which the kernel includes all fourth order contributions, but no higher order ones and which does not diverge in the stationary (zero-frequency) limit. This master equation can be derived using either the BR or RT (or NZ) approaches and there is no need to distinguish between these in the following discussion.

After this comparison on a formal level, we will continue in Sec. III with a comparison on the level of the individual contributions to the time-evolution kernel.

## Iii Diagrammatic representation and mapping between BR and RT

We now address the task of calculating all elements of the time evolution kernel in the stationary GME, Eq. (17). For our purposes, it will turn out to be advantageous to first work in the time-domain, i.e., to calculate , which we decompose into contributions of successive non-vanishing even orders in the tunnel coupling:

The section has a twofold aim. We first introduce the diagrammatic representation for the time evolution kernel and show how each BR contribution, obtained from Eq. (31), translates into a corresponding diagram in the RT approach based on Eq. (38). Apart from being of technical interest, this is of key importance for the discussion in Sec. IV, where the fourth order kernel elements incorporating corrections to the diagonal elements of the density matrix due to the non-diagonal elements are introduced. Secondly, we discuss the time-dependent part of the kernel contributions in the Schrödinger picture and its zero-frequency Laplace transform, on which the simplified calculation of the effective fourth order kernel in Sec. V relies.

In the conventional RT approach one starts by considering super matrix elements [cf. Eq. (18)] of the kernel ^{1}^{1}1See, however, the recently introduced superoperator formulation of the real-time approach Schoeller and Reininghaus (2009); Leijnse and Wegewijs (2008) which does not refer to matrix elements.
and one introduces a diagrammatic representation for the order parts of the kernel:

The diagram represents operators which act from the left and right on the dot density operator ,
inducing an irreducible time-evolution of a pair of initial states [associated with ] to corresponding final states [belonging to ].
Time thus increases from right to left.
The diagram can be considered as a Keldysh contour, i.e., running from , then continuing backward from , as indicated by the directed line on the upper and lower part.
On the upper (lower) part of the contour the time-ordering agrees with (is opposite to) the contour direction, indicating that operators acting from the right on the density operator come in inverted order concerning time.
This distinction is important for the diagram rules.
The shaded area indicates the sum of all contributions involving the product of tunnel operators , starting at time and ending at time (this thus yields a product of broadening elements, which we indicate with ).
Starting from the RT expression for the time evolution kernel, Eq. (38), simple rules are derived, given in App. B.2,
from which one can directly read off the analytical expression for the zero frequency Laplace transform of each diagram. Schoeller (1997)

Hence, the diagram contains not only the information about the contribution to the kernel in the interaction picture, but also to its Laplace transform . To make a distinction, we will use the convention that we mean contribution to whenever we explicitly write down time labels in the diagram. Otherwise it stands for the Laplace transformed expression, i.e. the contribution to (this is the case everywhere in the following except for Fig. 2).

In the BR approach one has to expand the superoperator expression Eq. (31) in commutators with dot and electrode operators. One then applies Wick theorem to integrate out the electrodes, resulting in cancellations of terms. Finally super-matrix elements are taken and the remaining expressions correspond term by term to the RT diagrams. To emphasize the close connection between the two approaches we now illustrate this explicitly by calculating second and fourth order contributions to the time evolution kernel in the BR approach. To this end, we first split the tunneling Hamiltonian (3) into two parts,

describing tunneling into () and out of () the dot:

(42a) | |||

(42b) |

Here the index numbers the time argument, with and corresponding to and , respectively. The summations over are implicitly understood. We insert this into Eq. (31), denoting :

(43a) | ||||

(43b) |

We next expand the multiple commutators and collect the fermionic operators of the leads, using that they anti-commute with the quantum dot operators, , where .

Using the cyclic property of the trace and Wick’s theorem, the average of the lead fermionic operators is expressed as a sum over products of pair-contractions. This results for Eq. (43a) and Eq. (43b) in the expressions listed on the left sides of Fig. 3. For the fourth order, the reducible contractions emerging from the last two terms in Eq. (31) cancel each other. Using the time-ordering of Fig. 2, Fig. 3 further gives on the right sides the respective RT diagram for each expression. The Hermitian conjugated terms, which have been omitted in the figures, correspond to diagrams which are vertically mirrored, i.e. all vertices on the upper contour have to be moved to the lower one and vice versa.

The translation from the BR expressions in the interaction picture into the RT diagrams works as follows. For each operator standing on the left (right) of , draw a vertex on the upper (lower) contour at time . For each contraction , requiring in order not to vanish, draw a contraction line between the vertices representing and . Notice that the ordering of the -operators in each contraction is consistently incorporated in the diagram: the pair of operators in the contraction have the same time-ordering as the corresponding operators, unless the earlier vertex of the two lies on the lower part of the contour (this follows from the cyclic permutation under the trace). Similarly, the sign of the operator expression, , is automatically contained in the diagram through the number of contractions (=order of perturbation theory), the number of crossing contraction lines , and the number of vertices on the lower contour .

The diagrams listed in Fig. 3 represent expressions which are summed over the indices . Terms with specific values of the ’s are represented by diagrams where the contraction lines are directed by an arrow, pointing towards the vertex corresponding to . Figure 4 shows an example for the third 4th order diagram in Fig. 3. From the diagrams it is explicitly clear that all contributions which where not canceled are irreducible: between the first and the last vertex at times , respectively , there is no time point at which the diagram could be separated into two parts without cutting a contraction line.

To obtain the stationary kernel Eq. (31) required in Eq. (15) we first need to transform back to the Schrödinger picture [cf. Eq. (9)] by inserting

(44a) | |||||

(44c) | |||||

For the further developments in this paper only the time-dependent part of the resulting expression, and its Laplace transform, is of importance. We can only factor out this part after taking super-matrix elements [cf. Eq. (18)] of the kernel contributions with respect to the energy-eigenstates of the quantum dot and insert complete sets of these states between all the vertex operators . The resulting expression is represented by a diagram labeled with these dot eigenstates on the contour, as illustrated in Fig. 4. We calculate its Laplace transform with respect to (), collecting for each time all energy contributions into one exponential with argument :

(45) |

Here the additional contains both the Laplace variable and the energy difference of the initial states on the upper and lower part of the contour. Transforming variables to the time-differences between vertices decouples the integrals, showing that the energies fully determine the time-evolution factor and its zero-frequency transform:

(46) |

This is the form of the zero frequency Laplace transform of the time-dependent factor only, as obtained from the diagram rules in App. B.2, which is the most convenient starting point for our analysis. The energy is obtained from a diagram by making a vertical cut through the diagram between times and , and then adding to / subtracting from the energies of all directed lines which hit this cut from the left / right. This includes the energies associated with contraction lines as well as the upper and lower part of the contour. Figure 4 demonstrates how this simple rule works for a specific fourth order diagram.

A crucial point for the rest of the paper is that diagrams which differ by breaking time-ordering between the contours, but keeping time-ordering within each contour, only differ by the arguments of the time-dependent exponentials in Eq. (46), i.e., in the . As App. B.2 shows, the products of TMEs and electrode distribution functions and the overall phase factor are identical. The simplifications discussed below are thus independent of these factors and their precise form needs no further discussion.

## Iv Coherences and non-secular corrections

For many simple quantum dot models, selection rules deriving from symmetries prevent the occupations of the dot states from coupling to the coherences. Whenever two states and of the system differ by some quantum number which is conserved in the total system (i.e. including the reservoirs), then their coherence does not enter into the calculation of the occupancies since for all states. The simplest example for such a quantum number is the electron number which is conserved for a quantum dot coupled to non-superconducting electrodes. The total spin-projection is also conserved for unpolarized or collinearly polarized electrodes. For non-collinear polarizations, however, inclusion of the coherences is crucial in order to capture spin-precession effects Braun et al. (2004, 2006). In a similar way orbital degeneracies have been shown to affect the occupations through the coherences. Begemann et al. (2008); Donarini et al. (2009) So in general, coherences cannot be neglected. Making no specific assumptions about the coherences, the only selection rule we enforce here is the one due to the conservation of total charge. In the above mentioned works, the tunneling is treated to lowest order and only non-diagonal elements between degenerate states are kept. The latter so-called secular approximation is usually phrased as neglecting the rapidly oscillating terms corresponding to coherences between non-degenerate states. Blum (1996)

In fourth order, however, an elimination of coherences between non-degenerate states in the stationary limit requires an expansion of the effective kernel for the occupations. Such an expansion is consistent in the sense that the derived effective kernel includes all contributions up to fourth order (while a simple neglect of non-secular coherences introduces serious errors of the order ). Leijnse and Wegewijs (2008); Koller (2010)

We start by decomposing the density matrix into a secular (energy diagonal) part and a non-secular (energy non-diagonal) part . Here contains all matrix elements between states with and all other elements (including the diagonal components, , corresponding to the populations). The cutoff should be chosen large compared to the tunnel broadening of the states, , the precise requirement being that it should be large enough that the next-to-leading order term in the expansion of Eq. (47) below is comparable to a sixth order term, and can thus be neglected. Our aim is to eliminate the non-secular coherences and include their effect as a correction to the kernel determining the secular part. To this end we write Eq. (15) in block matrix form,

where, see Eq. (17), the free evolution of the system involving is zero in the and blocks by definition. Solving for one obtains

(47) |

which obviously contains all orders in due to the inverse. Since by construction we can expand Eq. (47) in , finding that the lowest order term gives corrections to the secular density matrix of order and is thus all that should be kept in a consistent fourth order expansion. Inserting this in the equation for the secular part of the density matrix we obtain an effective stationary kinetic equation for the secular density matrix:

(48) |

with the effective fourth order kernel

(49) |

Here

(50) |

is the correction to the secular density matrix due to coherences between non-secular states. This makes explicit that when going beyond lowest order, the secular approximation is no longer valid and also coherences between non-secular states have to be accounted for. This was shown in Leijnse and Wegewijs (2008) for the special case where the secular part is diagonal. Here we extended the derivation to an arbitrary excitation spectrum, where the secular part may be non-diagonal, i.e., the effective equation is not a master-equation for occupancies. It should be noted that in this case the kernel must be calculated including the elements which couple to secular coherences.

For the developments of the present paper it is useful to introduce a diagrammatic representation of the non-secular correction . We first note that the inverse of is related to a diagram without tunneling lines by the diagram rules, see App. B.2, evaluated at zero frequency ():

Note that this ”free evolution” term is always finite since the expansion is only carried out in the non-secular subspace where and can always be dropped. Representing a general second order contribution diagrammatically as

the correction term due to coherences between non-secular states is given by