Counting statistics of cotunneling electrons

# Counting statistics of cotunneling electrons

Clive Emary Institut für Theoretische Physik, Technische Universität Berlin, D-10623 Berlin, Germany
July 30, 2019
###### Abstract

We describe a method for calculating the counting statistics of electronic transport through nanoscale devices with both sequential and cotunneling contributions. The method is based upon a perturbative expansion of the von Neumann equation in Liouvillian space, with current cumulants calculated from the resulting nonMarkovian master equation without further approximation. As application, we consider transport through a single quantum dot and discuss the effects of cotunneling on noise and skewness, as well as the properties of various approximation schemes.

###### pacs:
73.23.Hk, 73.23.-b, 73.63.Kv, 42.50.Lc

Cotunneling, the transfer of electrons via intermediate “virtual” states, can be an important mechanism in the transport of electrons through quantum dots (QDs) ave90 (); ave92 (). In the Coulomb blockade regime, sequential tunnelling processes are exponentially suppressed and, since it only suffers an algebraic suppression, cotunneling becomes the dominant current-carrying mechanism. Experimental interest in cotunneling has remained high from the earliest experiments on metallic grains gee90 () and large quantum dots pas93 (), through to more modern experiments on few-electron single- fra01 (); zum04 (); sch05 () and double- sig06 (); gus08 () quantum dots.

From a theoretical perspective, cotunneling refers to processes fourth-order in the coupling between the system and the leads. Such processes can be taken into account in a number of different ways, e.g. ave92 (); suk01 (); gol04 (); ped05 (); sch94 (); kon96 (); thi03 (); thi05 (). Most relevant here is the real-time diagrammatic approachsch94 (); kon96 (); thi03 (); thi05 (), in which higher-order tunnelling processes are incorporated into a master equation in a systematic fashion . This theory has been extensively developed and successfully applied to numerous transport problems: not just single QDs, but also double dots wey08DQD (), quantum dot spin-valves bra04 (), carbon nanotubes wey08CNT (), and QD-interferometers urb08 (). Whilst such higher-order calculations have typically been restricted to the stationary current, and more recently, the shotnoise thi03 (); thi05 (), much interest presently surrounds the full counting statistics (FCS) of the current, i.e. in current correlations beyond the second-order shotnoise lev93 (); yul99 (); bag03 (). The last few years has seen the advent of experiments capable of detecting the passage of single electrons through QD systems fuj06 (); gus06 (); Sukh07 (); fri07 (); fli09 () and the experimental determination of FCS. Recently, measurements of the 15th cumulant were reported for a single quantum dotfli09 ().

In this article we bring together several strands to investigate the influence of cotunneling on FCS. We derive a fourth-order master equation for the reduced density matrix of an arbitrary mesoscopic system using the Liouvillian-space perturbation theory of Ref. lei08 () and show how counting fields may be added in this formalism. Given that the resulting master equation is nonMarkovian, we employ the nonMarkovian formalism of Flindt et al. fli08 () to obtain expressions for the current cumulants.

The calculation of the FCS for a nonMarkovian ME with cotunneling was described in Ref. bra06 (), but the approach followed here is somewhat different. Braggio et al. bra06 () calculate the cumulant generating function (CGF), and thus all the cumulants, rigorously to fourth order in the dot-lead coupling. Here, we calculate the Liouvillian rigorously to fourth-order and make no further approximations in obtaining the cumulants. It is one of the aims of this paper to investigate the differences between the predictions of these two approaches.

We use the transport through a single QD as our example system. We study first the single resonant level (SRL) model. Exact solutions exist for this model and this allows an evaluation of various approximation schemes. We also study the effects of interaction on transport with an Anderson model. Of the higher-order cumulants, we focus here on the skewness as the first correlator beyond the shotnoise. We compare, both on a formal and a numerical level, with the work of Braggio et al bra06 (), and with the shotnoise results of Thielmann et al thi05 ().

## I Transport Model

We begin by specifying the general transport set-up under consideration here. The total Hamiltonian is composed of reservoir, system, and interaction parts. We write the system part in its diagonal basis , where is a many-body system state of electrons. We consider a set of reservoirs, labelled with an index that includes spin and any other relevant quantum numbers. We assume noninteracting reservoirs with Hamiltonian

 Hres=∑k,α(ωkα+μα)a†kαakα, (1)

where is the energy of the th mode in lead , is a lead annihilation operator, and we have included the chemical potential of lead , , at this point for convenience.

To ease book-keeping, we introduce a compact single index “” to denote the triple of indices . The first index indicates whether a reservoir operator is a creation or annihilation operator:

 a1=aξ1k1α1={a†k1α1,ξ1=+ak1α1,ξ1=−. (2)

Leaving sums implicit, the reservoir Hamiltonian reads

 Hres = (ωkα+μα)a+kαa−kα (3) = (ω1+μ1)a1a¯1δξ1,+.

In equilibrium, the reservoir electrons are distributed according to the Fermi function

 f(ω)=1eω/kBT+1, (4)

which, since we include the chemical potential in Eq. (1) and assume a uniform temperature, is the same for all reservoirs.

Single-electron tunnelling between system and reservoirs is described by the Hamiltonian

 V=∑kαmtkαma†kαdm+t∗kαmd†makα, (5)

where is the annihilation operator for single-particle level in the system, and is a tunnelling amplitude. We write this interaction as

 V = ξ1t1ma1jξ1m, (6)

with coefficients and , and system operators in the many-body system basis

 j+m=∑aa′⟨a|dm|a′⟩δ(Na−Na′+1)|a⟩⟨a′|, (7)

and . We have made explicit here the change in system charge induced by the operator. Although these operators only depend on , we label them with the full “” index for convenience: .

At time we posit a separable total density matrix:

 ρ(t=0)=ρS(0)ρeqres, (8)

with the system in arbitrary state and reservoirs in thermal equilibrium.

## Ii Liouville-Laplace space

We now construct the elements required to perform our perturbation calculation in Liouville-Laplace space. In this section and the next, we follow Refs.lei08 (); kor07 (); sch08 (). The total density matrix evolves according to the von Neumann equation:

 ˙ρ(t)=−i[H,ρ(t)]=Lρ(t), (9)

which defines the Liouvillian super-operator . This Liouvillian consists of three parts:

 L=Lres+LS+LV (10)

with , , and . We write the interaction Liouvillian as

 LV=−iξ1t1m∑pAp1Jp1m, (11)

where is a Keldysh index corresponding to the two parts of the commutator. Superoperators and are defined through their actions on arbitrary operator : For the reservoir, we have

 Ap1O={a1O,p=+Oa1,p=−, (12)

and analogously for the system

 Jp1mO={j1mO,p=+Oj1m,p=−. (13)

By organising the elements of density matrices into vectors, superoperators, such as the Liouvillian, take the form of matrices. This is a particularly convenient representation for the system Liouvillian. We write a general system density matrix, , as the vector , where the single index corresponds to the double such that the “ket” corresponds to . The action of the free system Liouvillian on vector is

 LS|ϕa⟩⟩≡−i[HS,|a1⟩⟨a2| ]=−iΔa|ϕa⟩⟩, (14)

where defines the Bohr frequencies. The vectors are therefore the right eigenvectors of Liouvillian . The left eigenvectors, , fulfil

 ⟨⟨a|LS=−iΔa⟨⟨ϕa| (15)

and together with the right eigenvectors form a bi-orthonormal set: . We have the completeness relation in Liouvillian space

 \mathbbm1=∑a|ϕa⟩⟩⟨⟨ϕa|. (16)

In general, it is important to make the distinction between left and right eigenvectors because an arbitrary super-operator, in particular the effective system Liouvillian, will not be Hermitian, and the left and right eigenvectors are therefore not adjoint.

## Iii Effective Liouvillian

With the definition of the Laplace transform

 ρ(z)≡∫∞0dte−ztρ(t), (17)

equation (9) yields the solution

 ρ(z)=1z−Lρ(0). (18)

Tracing Eq. (18) over reservoir degrees of freedom, results in an expression for the reduced density matrix of the system that we write

 ρS(z)=1z−W(z)ρS(t0) (19)

where is the nonMarkovian effective dot Liouvillian. This we write as

 W(z)=LS+Σ(z), (20)

with describing the free evolution of the system and the self-energy or “memory kernel” arising from coupling with the leads.

In the perturbative approach pursued here, the memory kernel is calculated as the series where corresponds to the number of interaction Liouvillians incorporated in that term. Tunnelling is governed by the rates

 Γm2m1ξ1α1(ω)≡2π∑k1t¯1m2t1m1δ(ω−ωkα), (21)

the diagonal elements of which are the familiar Fermi golden rule rates

 Γm1m1ξ1α1(ω)≡2π∑k1|t1m1|2δ(ω−ωkα). (22)

In these terms, the expansion of is seen as an expansion in small parameter , such that is order . In the current work, we expand up to fourth order in the coupling Hamiltonian (second order in ), such that

 Σ(z)≈Σ(2)(z)+Σ(4)(z). (23)

The first term describes sequential tunnelling and the second cotunneling.

Details of the calculation of the memory kernel terms are given in Appendix A. Assuming a constant tunnelling density of states , the sequential term reads

 Σ(2)(z) = Jp22m2−p2z−iξ2(ω2+μ2)−LSJp11m1 (24) ×t2m2t1m1γp2p121,

where is an equilibrium reservoir correlation function which evaluates as

 γp2p121=δ2¯1p1f(−ξ1p1ω1). (25)

We then obtain

 Σ(2)(z) = −p1p2Jp2¯1m2|ϕa⟩⟩⟨⟨ϕa|Jp11m1 (26) ×Γm2m1ξ1α1I(2)a(z;ξ1,p1,μ1),

with the regularised integral

 I(2)a(z=O+−iϵ) = 12f(p1(Δa+ξ1μα1−ϵ)) +ip12πϕ(p1(Δa+ξ1μα1−ϵ)),

in terms of the function

 ϕ(λ) = ReΨ(12+iλ2πkBT)−lnD2πkBT (27)

with the digamma function.

The cotunneling term has two contributions: “direct” and “exchange”, such that . The direct part is given by

 Σ(4D)(z) = p4p1Jp4¯1m4|ϕa⟩⟩⟨⟨ϕa|Jp3¯2m3|ϕa′⟩⟩ (28) ×⟨⟨ϕa′|Jp22m2|ϕa′′⟩⟩⟨⟨ϕa′′|Jp11m1 ×Γm4m11Γm3m22 ×IDaa′a′′(z;ξ1,ξ2,p1,p2,μ1,μ2),

and the exchange

 Σ(4X)(z) = −p4p1Jp4¯2m4|ϕa⟩⟩⟨⟨ϕa|Jp3¯1m3|ϕa′⟩⟩ (29) ×⟨⟨ϕa′|Jp22m2|ϕa′′⟩⟩⟨⟨ϕa′′|Jp11m1 ×Γm3m11Γm4m22 ×IXaa′a′′(z;ξ1,ξ2,p1,p2,μ1,μ2).

These fourth-order integrals, and , are discussed in the Appendix.

## Iv Full Counting Statistics

The density matrix of Eq. (19) is the Laplace-transform of the solution to the nonMarkovian master equation FN_nonMark ()

 ˙ρS(t)=∫t0dt′W(t−t′)ρ(t′). (30)

In contrast to e.g. Ref. lei08 (), we make no further approximations at this point — in particular, we do not make Markov approximation — but rather seek to calculate the FCS of Eq. (30) as it stands. This requires that we introduce counting fields in the appropriate places in the Liouvillian, a process which yields the “-resolved” Liouvillian bag03 ().

In the current Liouville scheme, this can be achieved by replacing each contraction in the memory kernel by the counting-field-dependent analogue , which we define as

 γp2p121(χ)=γp2p121exp[isα1ξ1(p1−p22)χα1], (31)

with a factor given by the sign-convention for current flow in lead . In the following we will only count in a single lead, for which we choose . To see that adds a counting field at the correct points, consider, for example, the case with . Then, for , the contraction is proportional to trace of , a state with one more electron in lead than itself. On the other hand, for , we require the trace of , a state with the same number of electrons as . The required counting-field factors are and , respectively, as given by Eq. (31).

Once in possession of the “-resolved” Liouvillian, the cumulant generating function (CGF) is obtained from the solution of the equation

 z0−λ0(χ;z0)=0, (32)

where is the eigenvalue of that develops adiabatically from zero as is increased from zero jau05 (). In the Markovian case, is independent of , and the CGF is simply . The nonMarkovian case is less straightforward, however. We follow here the approach of Ref. fli08 (), which uses Eq. (32) to derive expressions for the cumulants themselves, by-passing an explicit evaluation of the CGF itself. This approach can deliver the cumulants up to, in principle, arbitrary order ( in Ref. fli04 ()) and without any further approximation.

We first define

 J(χ,ϵ)=W(χ,z=0+−iϵ)−W(χ=0,z=0+), (33)

along with the derivatives

 J′=∂χJ∣∣χ,ϵ→0;˙J=∂ϵJ|χ,ϵ→0, (34)

and analogously for higher-orders. We define the left and right null-vectors of via

 W(0,0+)|ψ0⟩⟩=⟨⟨ψ0|W(0,0+)=0, (35)

which we assume to be unique. The vector corresponds to the stationary density matrix of the system, and multiplication with , corresponds to taking the trace over system statesjau05 (). We define the stationary state “expectation value” , and the projectors and . Finally, we require the pseudo-inverse

 R(ϵ)=Q1iϵ+L(0,0+−iϵ)Q. (36)

From Refs. fli08 (); flindtthesis (), the first three current cumulants are

 ⟨⟨I1⟩⟩ = ⟨⟨I1⟩⟩m (37) ⟨⟨I2⟩⟩ = ⟨⟨I2⟩⟩m+2⟨⟨I1⟩⟩m⟨⟨˙J′−J′R˙J⟩⟩ (38) ⟨⟨I3⟩⟩ = ⟨⟨I3⟩⟩m−3⟨⟨I2⟩⟩2⟨⟨I1⟩⟩m(⟨⟨I2⟩⟩m−⟨⟨I2⟩⟩) (39) −3i⟨⟨I1⟩⟩m⟨⟨˙J′′−2˙J′RJ′−J′′R˙J⟩⟩ +6i⟨⟨I1⟩⟩m⟨⟨J′R(R˙JPJ′+˙J′)⟩⟩ −6i⟨⟨I1⟩⟩m⟨⟨J′R(J′R˙J+˙JRJ′)⟩⟩ +3i(⟨⟨I1⟩⟩m)2⟨⟨¨J′+2J′R˙JR˙J⟩⟩ −3i(⟨⟨I1⟩⟩m)2⟨⟨J′R¨J+2˙J′R˙J⟩⟩,

where

 i⟨⟨I1⟩⟩m = ⟨⟨J′⟩⟩ (40) i2⟨⟨I2⟩⟩m = ⟨⟨J′′−2J′RJ′⟩⟩ (41) i3⟨⟨I3⟩⟩m = ⟨⟨J′′′−3J′′RJ′−3J′RJ′′⟩⟩ (42) −6⟨⟨J′R(RJ′P−J′R)J′⟩⟩,

are the cumulants in the Markov approximation. In these expressions, it is understood that the pseudo-inverse is evaluated at . Although it is only practicable to explicitly write down the cumulants up to third order, the high-order cumulants can be obtained recursively fli08 ().

Reference bra06 () took a different approach to calculating the cumulants. There it was assumed that is known to a given order in some small parameter, and the CGF is then calculated to the same order. For problems such as that considered here, this means that the CGF, and hence all the cumulants, are calculated rigorously up to order . This is to be contrasted with the above cumulants which, if expanded, have contributions at all orders in . Whilst the method of Ref. bra06 () may seem more consistent, there are two good reasons why the approach described here might be preferable. Firstly, from Ref. gur96 () we know that in the infinite bias limit, the effective Liouvillian is given exactly by plus the rate part of . In order to recover the FCS correctly in this limit then, no further approximations should be made when calculating the cumulants FNDQD (). Secondly, Ref. lei08 () makes the point that, in certain circumstances, by treating incoming and outgoing processes unequally, a strict order-by-order approach can lead to unphysical results for the current. In the next section we shall compare these two methods directly.

Before doing so, we note that comparing the foregoing expressions for the current and the shotnoise with those of, e.g., Ref. thi05 (), we can identify the derivatives of with respect to and with the various “current super-operator blocks” of the real-time diagrammatic approach. For example, differentiating once with respect to and setting yields a super-operator like , but with an additional forefactor . Under the summation, the contribution cancels, leaving a forefactor . This forefactor is the same as arises in replacing a single tunnel vertex with a current vertex in the self-energy, which is the recipe for obtaining the current super-operator block used to calculate the stationary current in the diagrammatic approach. Using a very different approach, we thus reproduce the real-time current and shotnoise expressions. The advantage of the present method is that it is now easy to obtain to the higher cumulants, whereas with the diagrammatic approach, this requires some effort.

## V Transport through a single quantum dot

We model the transport through a single Zeeman-split level with the Anderson Hamiltonian and61 ()

 H=∑σϵσd†σdσ+Un↑n↓+Hres+V, (43)

where is the energy of a spin- electron in the dot and is the interaction energy. The reservoir and interaction Hamiltonians are as Eq. (1) and Eq. (5), with index including both lead () and spin index. In the limit of large level-splitting we can address one and only one Zeeman level. We then recover the SRL model, the transport properties of which can be obtained exactly from scattering theory but92 (); bla00 (). The SRL thus provides a useful benchmark against which to compare our approximate methods.

We will discuss results obtained in several different approximate schemes. We denote as “O(4)” the results obtained by calculating the cumulants as outlined in the previous section with the fourth-order effective Liouvillian containing both sequential and cotunneling terms. The second-order “O(2)” solution is obtained in the same way, but with sequential terms only. We also consider a scheme in which we expand the cumulants and truncate at fourth order. In this way we recover the FCS results of Braggio et al bra06 () and the shotnoise results of Thielmann et al thi05 (). This approach we label as “O(4) trunc”. Finally, we compare with results in the Markovian approximation, which we label with “Mark.”.

We calculated results with and without the level-renormalisation parts of the fourth-order self-energy (integrals and , and and of the Appendix). Their contribution was found to be negligible in all cases studied here. In the results presented below, these parts of the self-energy have been neglected. This reduces considerably the computational effort involved since the double-principal-part integrals ( and ) must be evaluated numerically.

### v.1 Single resonant level

The calculation of the first three cumulants in the scattering approach is discussed in Appendix B. In the infinite bias limit, we have jon96 ()

 ⟨I⟩ = ΓLΓRΓ;S=⟨I⟩Γ2L+Γ2RΓ2; S(3) = ⟨I⟩Γ4(Γ4L−2Γ3LΓR+6Γ2LΓ2R−2ΓLΓ3R+Γ4R),

or , and for symmetric coupling, .

Figure 1 shows the current and shotnoise as a function of applied bias for different values of the coupling . A step occurs at ( here) as the level enters the transport window. For agreement between the our O(4) fourth-order calculation and the exact result is excellent across the whole bias range. For greater couplings, , deviation from the exact solution is seen in the shotnoise around the top of the step which signals the start of the break down of our approach. The difference between the second-order Markovian and the exact solution is stark. In the CB regime ( here) the sequential current is almost totally suppressed, but the cotunneling current is still considerable. The O(2) Mark. solution also shows significant error in the high-bias () regime, which arises largely from the Markovian approximation.

Figure 2 plots the skewness which show a pronounced undulation at the onset. The O(2) Markovian solution provides only the coarsest description of this behaviour, whereas it is reproduced by the O(4) solution. For , the quantitative agreement with the exact result is good. Nevertheless, for a given coupling, the error is larger for skewness than for the noise. This is not too surprising, since we expect higher-order correlators to be sensitive to higher-order processes, which are of course neglected here. The implication is that for a given coupling, only cumulants up to a certain order can reliably be calculated from any approximate effective Liouvillian FNinf (). Figures 3 and 4 compare the O(4) and O(4) trunc results at . We choose this value to highlight the differences between the two solutions which, for smaller couplings, are negligible.

Near the top of the step (Fig. 3a and Fig. 4a), the O(4) and O(4) trunc solutions are noticeably different, with our O(4) results significantly closer to the exact result. Deep in the CB regime (Fig. 3b), the two approximate solutions differ once more, but this time the trunc solution is more accurate. The difference is small (a difference in of near zero bias); it plays, however, a disproportionate role in determining the Fano factors in the CB regime as Fig. 3c and Fig. 4b show. From the exact solution, we know that noise Fano factor diverges at low bias (fluctuation dissipation theorem), whereas the skewness Fano factor tends to unity. In this regime, both Fano factors are reproduced better by trunc solution than by solution O(4). It should be borne in mind that, in obtaining the Fano factors in this region, one is dividing one very small quantity by another and thus even small absolute errors can effect Fano factors quite dramatically. As a warning not to take the finer details of the Fano factor too seriously, we observe that the O(2) Mark. solution in the CB regime reproduces the exact Fano factors better than any of the fourth-order results. This is deceptive since the individual current and shotnoise obtained with this method are vastly different from the exact results.

From this comparison we obtain a degree of confidence in the cumulants of up to at least third order for . Our O(4) solution appears to perform better in the region where levels are crossing the chemical potentials, whereas the O(4) trunc solution gives better results in the CB regime. In this regime, the transport properties are effectively Markovian, which can be demonstrated by comparing fourth-order solutions with and without the Markov approximation..

### v.2 Anderson Model

The current, shotnoise and Fano factor of the Anderson model with cotunneling were investigated in Ref. thi05 (), and, as Fig. 5 shows, the present calculation broadly reproduces these results. The situation in which the lower dot level lies below the transport window is of particular interest. As observed in Ref. thi05 (), increasing the applied bias results in a large peak in the Fano factor around the point where the upper dot level enters the transport window. The peak exists in the sequential tunnel limit, but its height, width, and location are markedly altered by inelastic cotunneling processes. No peak occurs in the shotnoise itself; only in the Fano factor is this feature visible. Figure 6 shows our results for the skewness in this situation. It is immediately clear that the skewness Fano factor also shows a peak, and that this is even more pronounced than that of the noise. Furthermore, the skewness itself exhibits a sharp peak, as inset Fig. 6b shows. As for the noise Fano factor, the presence of cotunneling significantly reduces the height and overall area of the peak.

This superPoissonian behaviour indicates a significantly bunched electron flow, which can nicely be explained with the dynamical channel blockade model of bel04 (), in which a single level (here, the lower) is but weakly coupled to the collector. In the simple sequential picture of Ref. bel04 (), the shotnoise and skewness Fano factors are predicted to be and , where is parameter corresponding to the number of ways in which the dot can be filled. With (corresponding to three ways of filling the dot: from the left and right into the lower level, and from the left only into the upper level), we obtain and , which are almost exactly the values obtained by our sequential O(2) results at the tops of the peaks. Cotunneling reduces the heights of peaks, and good agreement with the dynamical channel blockade model can be obtained with the choice .

Both noise and skewness figures show results for O(4) and O(4) trunc solutions. At high bias, these solutions agree closely with one another and both predict the same position and widths of the Fano factor peaks. However, the heights of the peaks are given differently in the two approaches, with the O(4) trunc peaks being somewhat higher. Differences between the two solutions are most pronounced at low bias, when the Coulomb blockade is in effect. For example, the O(4) solution predicts a small subPoissonan dip before the superPoissonian peak, which is absent in the O(4) trunc results. From our studies of the SRL model, we expect the O(4) trunc prediction to be more accurate in this regime, and this is borne out by the fact that the skewness Fano factor in the O(4) calculation does not tend to unity with decreasing bias. Conversely, based on the SRL results, we expect the height of the Fano factor peaks to be better described by the O(4) solution, i.e. the lower of the two values.

## Vi Conclusions

We have described a method for calculating the counting statistics of an arbitrary mesoscopic system that takes into account both sequential tunnelling and cotunneling of electrons. This is achieved by performing a perturbative expansion of the von Neumann equation in Liouville-Laplace space and adding counting fields to the reservoir contractions. Current cumulants are then obtained without further approximation with the pseudo-inverse approach. In principle, cumulants of arbitrary order can be calculated in this fashion. However, as we expect higher-order cumulants to be sensitive to higher-order tunnel processes, there will be an inevitable reduction in accuracy as the order of the cumulant increases. Use of the pseudo-inverse means that this method is applicable to systems of large size, unlike methods which explicitly require the eigenvalue of the effective Liouvillian.

We have studied here transport through a single quantum dot with both sequential and cotunneling contributions. Of particular interest is the single-resonant level model as it provides an exact case against which we can compare. Good agreement with the exact results was found for both shotnoise and skewness up to ratio between coupling rate and reservoir temperature of . Moreover, we also compared our results with those obtained from a rigorous expansion of the cumulants up to fourth order in the tunnel coupling. From this comparison we conclude that for intermediate biases, and in particular when system levels lie near the chemical potentials of the leads, the approach described here gives better results. In the Coulomb blockade regime, however, the rigorous fourth-order approach is better, but here the dynamics are effectively Markovian. In general, the difference between these two different fourth-order approximations is far less than that between fourth and second order approximations, and decreases with decreasing system-reservoir coupling.

Future work includes the study of transport models with internal quantum degrees of freedom, such as the double quantum dot, to understand how the interplay of cotunneling and internal dynamics effects counting statistics.

## Appendix A Derivation of effective system Liouvillian

Following Refs. lei08 (); kor07 (); sch08 (), we start by defining the system operators

 gkα = ∑mtkαmjm, (44)

such that the interaction Hamiltonian of Eq. (5) can be written . Correspondingly, in Liouville space we have

 LV=−iξ1∑ppσpAp1Gp1 (45)

with as before, and defined via

 Gp1O=σp{g1O,p=+−Og1,p=−. (46)

The object is a dot-space superoperator with matrix elements sch08 ()

 (σp)ss′,¯s¯s′=δs¯sδs′¯s′{1,Ns−Ns′= evenp,Ns−Ns′= odd, (47)

where, is the number of electrons in state . Note that .

The reduced density matrix of the dot is given by tracing out the electron reservoirs

 ρS(z)=TrR{ρ(z)}=TrR{1z−Lρ(0)}. (48)

This we expand in powers of to obtain

 ρS(z)=TrR{(Ω0(z)+Ω0(z)LVΩ0(z)+…)ρ(0)} (49)

with free propagator . With substitution of Eq. (45), a typical term of the expansion Eq. (49) reads

 (−i)n(n∏l=1ξlpl)TrR{Ω0(z)σpnApnnGpnn… …σp1Ap11Gp11Ω0(z)ρ(t0)}. (50)

Evaluating the action of the superoperators we obtain a factor and, as the operators also contain , they evaluate at different positions in the chain as

 GpllO = (pl)l+1{glOOgl. (51)

Our typical term then looks like

 (−i)n(all∏lξlpl)(odd∏l′pl′) ×TrR{Ω0(z)ApnnGpnn…Ap11Gp11Ω0(z)ρ(t0)}, (52)

and the next task is to separate dot and reservoir degrees of freedom. For this we can use the dot-reservoir superoperator commutation relation . We will also need the following reasons: , and with

 x1 = −iξ1(ω1+μα1). (53)

The commutation of the dot-operators through the free propagators therefore changes the argument of the propagator:

 Ap1Ω0(z)=Ω0(z+x1)Ap1. (54)

Bringing all the operators to the right of the operators generates a factor which exactly cancels with the first product in Eq. (52). Our term becomes

 (−i)n(odd∏l′pl′)TrR{ΩS(z)GpnnΩS(zn−1)Gpn−1n−1… …Gp22ΩS(z1)Gp11ΩS(z)Apnn…Ap11ρ(0)} (55)

with free dot propagator

 ΩS(z)=1z−LS, (56)

and .

The the reservoir expectation values, , can be evaluated with the rules of Wick’s theorem in Liouville space, which read lei08 ()

• decompose into pair-contractions

• add minus sign for each interchange of

• omit factor arising from super-operators

• each pair contraction then contributes a factor

 ⟨Ap22Ap11⟩=γp2p121=δ2¯1p1fα1(−ξ1p1ω1). (57)

With these rules, our typical term becomes

 (−i)nΩS(zn)GpnnΩS(zn−1)Gpn−1n−1… …Gp22ΩS(z1)Gp11ΩS(z)ρS(t0) ×(∑decomps(−1)NP∏γij), (58)

where the last factor indicates a sum over all pair decompositions with the relevant Wick sign .

Comparison of this expression with Eq. (19) and Eq. (20) allows us to identify the self-energy as

 Σ(z) = ×GpnnΩS(zn−1)Gpn−1n−1…Gp22ΩS(z1)Gp11,

where the sum is over irreducible contractions only.

### a.1 Second order

At second order, there is only one contraction, and we have

 Σ(2)(z) = Gp22−1z+x2−LSGp11γp2p121 = Gp2¯1−p1fα1(−ξ1p1ω1)z+iξ1(ω1+μα1)−LSGp11 = −Gp2¯1|ϕa⟩⟩⟨⟨ϕa|Gp11p1fα1(−ξ1p1ω1)z+iξ1(ω1+μα1)+iΔa.

Re-expressing the superoperators in terms of superoperators and tunnel amplitudes, we have

 Σ(2)(z) = −p1p2Jp2¯1m2|ϕa⟩⟩⟨⟨ϕa|Jp11m1 ×t¯1m2t1m1fα1(−ξ1p1ω1)z+iξ1(ω1+μα1)+iΔa.

With rates defined as in Eq. (21) with the constant tunnelling density of states approximation, Eq. (LABEL:SIG2tt) becomes

 Σ(2)(z) = −p1p2Jp2¯1m2|ϕa⟩⟩⟨⟨ϕa|Jp11m1 (60) ×Γm2m