Counting statistics of cotunneling electrons
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.LcCotunneling, 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 currentcarrying 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 fewelectron single fra01 (); zum04 (); sch05 () and double sig06 (); gus08 () quantum dots.
From a theoretical perspective, cotunneling refers to processes fourthorder 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 realtime diagrammatic approachsch94 (); kon96 (); thi03 (); thi05 (), in which higherorder 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 spinvalves bra04 (), carbon nanotubes wey08CNT (), and QDinterferometers urb08 (). Whilst such higherorder 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 secondorder 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 fourthorder master equation for the reduced density matrix of an arbitrary mesoscopic system using the Liouvillianspace 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 dotlead coupling. Here, we calculate the Liouvillian rigorously to fourthorder 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 higherorder 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 setup 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 manybody 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
(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 bookkeeping, 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:
(2) 
Leaving sums implicit, the reservoir Hamiltonian reads
(3)  
In equilibrium, the reservoir electrons are distributed according to the Fermi function
(4) 
which, since we include the chemical potential in Eq. (1) and assume a uniform temperature, is the same for all reservoirs.
Singleelectron tunnelling between system and reservoirs is described by the Hamiltonian
(5) 
where is the annihilation operator for singleparticle level in the system, and is a tunnelling amplitude. We write this interaction as
(6) 
with coefficients and , and system operators in the manybody system basis
(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:
(8) 
with the system in arbitrary state and reservoirs in thermal equilibrium.
Ii LiouvilleLaplace space
We now construct the elements required to perform our perturbation calculation in LiouvilleLaplace space. In this section and the next, we follow Refs.lei08 (); kor07 (); sch08 (). The total density matrix evolves according to the von Neumann equation:
(9) 
which defines the Liouvillian superoperator . This Liouvillian consists of three parts:
(10) 
with , , and . We write the interaction Liouvillian as
(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
(12) 
and analogously for the system
(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
(14) 
where defines the Bohr frequencies. The vectors are therefore the right eigenvectors of Liouvillian . The left eigenvectors, , fulfil
(15) 
and together with the right eigenvectors form a biorthonormal set: . We have the completeness relation in Liouvillian space
(16) 
In general, it is important to make the distinction between left and right eigenvectors because an arbitrary superoperator, 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
(17) 
equation (9) yields the solution
(18) 
Tracing Eq. (18) over reservoir degrees of freedom, results in an expression for the reduced density matrix of the system that we write
(19) 
where is the nonMarkovian effective dot Liouvillian. This we write as
(20) 
with describing the free evolution of the system and the selfenergy 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
(21) 
the diagonal elements of which are the familiar Fermi golden rule rates
(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
(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
(24)  
where is an equilibrium reservoir correlation function which evaluates as
(25) 
We then obtain
(26)  
with the regularised integral
in terms of the function
(27) 
with the digamma function.
The cotunneling term has two contributions: “direct” and “exchange”, such that . The direct part is given by
(28)  
and the exchange
(29)  
These fourthorder integrals, and , are discussed in the Appendix.
Iv Full Counting Statistics
The density matrix of Eq. (19) is the Laplacetransform of the solution to the nonMarkovian master equation FN_nonMark ()
(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 countingfielddependent analogue , which we define as
(31) 
with a factor given by the signconvention 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 countingfield 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
(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, bypassing 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
(33) 
along with the derivatives
(34) 
and analogously for higherorders. We define the left and right nullvectors of via
(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 pseudoinverse
(36) 
From Refs. fli08 (); flindtthesis (), the first three current cumulants are
(37)  
(38)  
(39)  
where
(40)  
(41)  
(42)  
are the cumulants in the Markov approximation. In these expressions, it is understood that the pseudoinverse is evaluated at . Although it is only practicable to explicitly write down the cumulants up to third order, the highorder 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 orderbyorder 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 superoperator blocks” of the realtime diagrammatic approach. For example, differentiating once with respect to and setting yields a superoperator 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 selfenergy, which is the recipe for obtaining the current superoperator block used to calculate the stationary current in the diagrammatic approach. Using a very different approach, we thus reproduce the realtime 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 Zeemansplit level with the Anderson Hamiltonian and61 ()
(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 levelsplitting 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 fourthorder effective Liouvillian containing both sequential and cotunneling terms. The secondorder “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 levelrenormalisation parts of the fourthorder selfenergy (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 selfenergy have been neglected. This reduces considerably the computational effort involved since the doubleprincipalpart 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 ()
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) fourthorder 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 secondorder 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 highbias () 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 higherorder correlators to be sensitive to higherorder 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 fourthorder 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 fourthorder 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 LiouvilleLaplace space and adding counting fields to the reservoir contractions. Current cumulants are then obtained without further approximation with the pseudoinverse approach. In principle, cumulants of arbitrary order can be calculated in this fashion. However, as we expect higherorder cumulants to be sensitive to higherorder tunnel processes, there will be an inevitable reduction in accuracy as the order of the cumulant increases. Use of the pseudoinverse 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 singleresonant 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 fourthorder approach is better, but here the dynamics are effectively Markovian. In general, the difference between these two different fourthorder approximations is far less than that between fourth and second order approximations, and decreases with decreasing systemreservoir 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
(44) 
such that the interaction Hamiltonian of Eq. (5) can be written . Correspondingly, in Liouville space we have
(45) 
with as before, and defined via
(46) 
The object is a dotspace superoperator with matrix elements sch08 ()
(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
(48) 
This we expand in powers of to obtain
(49) 
with free propagator . With substitution of Eq. (45), a typical term of the expansion Eq. (49) reads
(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
(51) 
Our typical term then looks like
(52) 
and the next task is to separate dot and reservoir degrees of freedom. For this we can use the dotreservoir superoperator commutation relation . We will also need the following reasons: , and with
(53) 
The commutation of the dotoperators through the free propagators therefore changes the argument of the propagator:
(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
(55) 
with free dot propagator
(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 paircontractions

add minus sign for each interchange of

omit factor arising from superoperators

each pair contraction then contributes a factor
(57)
With these rules, our typical term becomes
(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 selfenergy as
where the sum is over irreducible contractions only.
a.1 Second order
At second order, there is only one contraction, and we have
Reexpressing the superoperators in terms of superoperators and tunnel amplitudes, we have
With rates defined as in Eq. (21) with the constant tunnelling density of states approximation, Eq. (LABEL:SIG2tt) becomes
(60)  