Counting statistics of transport through Coulomb blockade nanostructures: High-order cumulants and non-Markovian effects

Counting statistics of transport through Coulomb blockade nanostructures:
High-order cumulants and non-Markovian effects

Christian Flindt Department of Physics, Harvard University, 17 Oxford Street, Cambridge, MA 02138, USA    Tomáš Novotný Department of Condensed Matter Physics, Faculty of Mathematics and Physics, Charles University, Ke Karlovu 5, 12116 Prague, Czech Republic    Alessandro Braggio CNR-SPIN, Dipartimento di Fisica, Università di Genova, Via Dodecaneso 33, 16146 Genova, Italy    Antti-Pekka Jauho Dept. of Micro og Nanotechnology, DTU Nanotech, Technical University of Denmark, Building 345east, 2800 Kongens Lyngby, Denmark Aalto University, Department of Applied Physics, P. O. Box 11100, FI-00076 AALTO, Finland
July 5, 2019

Recent experimental progress has made it possible to detect in real-time single electrons tunneling through Coulomb blockade nanostructures, thereby allowing for precise measurements of the statistical distribution of the number of transferred charges, the so-called full counting statistics. These experimental advances call for a solid theoretical platform for equally accurate calculations of distribution functions and their cumulants. Here we develop a general framework for calculating zero-frequency current cumulants of arbitrary orders for transport through nanostructures with strong Coulomb interactions. Our recursive method can treat systems with many states as well as non-Markovian dynamics. We illustrate our approach with three examples of current experimental relevance: bunching transport through a two-level quantum dot, transport through a nano-electromechanical system with dynamical Franck-Condon blockade, and transport through coherently coupled quantum dots embedded in a dissipative environment. We discuss properties of high-order cumulants as well as possible subtleties associated with non-Markovian dynamics.

02.50.Ey, 03.65.Yz, 72.70.+m, 73.23.Hk

I Introduction

Electron transport through nanoscale structures is a stochastic process due to the randomness of the individual tunneling events. Quantum correlations and electron-electron interactions can strongly influence the transport process and thus the statistics of transferred charges. Full counting statisticsLevitov and Lesovik (1993); Levitov et al. (1996); Nazarov (2003) concerns the distribution of the number of transferred charge, or equivalently, all corresponding cumulants (irreducible moments) of the distribution. Conventional transport measurements have focused on the first cumulant, the mean current, and in some cases also the second cumulant, the noise.Blanter and Büttiker (2000) Higher order cumulants, however, reveal additional information concerning a variety of physical phenomena, including quantum coherence, entanglement, disorder, and dissipation.Nazarov (2003) For example, non-zero higher-order cumulants reflect non-Gaussian behavior. Counting statistics in mesoscopic physics has been a subject of intensive theoretical interest for almost two decades, but recently it has also gained considerable experimental interest: in a series of experiments,Reulet et al. (2003); Bomze et al. (2005); Bylander et al. (2005); Fujisawa et al. (2006); Gustavsson et al. (2006); Fricke et al. (2007); Timofeev et al. (2007); Gershon et al. (2008); Gabelli and Reulet (2009); Flindt et al. (2009); Gustavsson et al. (2009); Fricke et al. (2010, 2010) high order cumulants and even the entire distribution function of transferred charge have been measured, clearly demonstrating that counting statistics now has become an important concept also in experimental physics.

The theory of counting statistics was first formulated by Levitov and Lesovik for non-interacting electrons using a scattering formalism.Levitov and Lesovik (1993); Levitov et al. (1996) Subsequent works have focused on the inclusion of interaction effects in the theory.Nazarov (1999); Kindermann and Nazarov (2003) In one approach, Coulomb interactions are incorporated via Markovian (generalized) master equations as originally developed by Bagrets and Nazarov.Bagrets and Nazarov (2003) This approach is often convenient when considering systems with strong interactions, e.g., Coulomb-blockade structures. More recent developments include theories for finite-frequency counting statistics,Emary et al. (2007) conditional counting statistics,Sukhorukov et al. (2007) connections to entanglement entropyKlich and Levitov (2009) and to fluctuation theorems,Förster and Büttiker (2008); Esposito et al. (2009) and extensions to systems with non-Markovian dynamics.Braggio et al. (2006); Flindt et al. (2008); Schaller et al. (2009) The last topic forms the central theme of this paper.

We have recently published a series of papers on counting statistics.Flindt et al. (2005a); Braggio et al. (2006); Flindt et al. (2008) Previous methods for evaluating the counting statistics of systems described by master equations had in practice been limited to systems with only a few states, and in Ref. Flindt et al., 2005a we thus developed techniques for calculating the first few cumulants of Markovian systems with many states, for example nano-electromechanical systems.Flindt et al. (2004) In Ref. Braggio et al., 2006, Braggio and co-workers generalized the approach by Bagrets and Nazarov by including non-Markovian effects that may arise for example when the coupling to the electronic leads is not weak. The methods presented in these papers were subsequently unified and extended in Ref. Flindt et al., 2008, where we presented a general approach to calculations of cumulants of arbitrary order for systems with many states as well as with non-Markovian dynamics. The aim of the present paper is to provide a detailed derivation and description of this method, which recently has been used in a number of works,Zedler et al. (2009); Emary (2009); Emary et al. (2009); Urban and König (2009); Lindebaum et al. (2009); Zhong and Cao (2009); Domínguez et al. (2010) as well as to illustrate its use with three examples of current experimental relevance.

The paper is organized as follows: In Sec. II we introduce the generic non-Markovian generalized master equation (GME) which is the starting point of this work. The GME describes the evolution of the reduced density matrix of the system, which has been resolved with respect to the number of transferred particles. Memory effects due to the coupling to the environment as well as initial system-environment correlations are included in the GME. Within this framework it is possible to calculate the finite-frequency current noise for non-Markovian GMEsFlindt et al. (2008) as we will discuss in future works. Section II concludes with details of the superoperator notation used throughout the paper.

In Sec. III we develop a theory for the zero-frequency cumulants of the current. The cumulant generating function (CGF) is determined by a single dominating pole of the resolvent of the memory kernel, and its derivatives with respect to the counting field evaluated at zero counting field yield the cumulants of the current. Even in the Markovian case it is difficult to determine analytically the dominating pole and in many cases one would have to find it numerically. Numerical differentiation, however, is notoriously unstable, and often one can only obtain accurate results for the first few derivatives with respect to the counting field, i. e. the cumulants. In order to circumvent this problem, we develop a numerically stable recursive scheme based on a perturbation expansion in the counting field. The scheme enables calculations of zero-frequency current cumulants of very high orders, also for non-Markovian systems. Some notes on the evaluation of the cumulants are presented, with the more technical numerical details deferred to App. A.

Section IV gives a discussion of the generic behavior of high-order cumulants. As some of us have recently shown,Flindt et al. (2009) the high-order cumulants for basically any system (with or without memory effects) are expected to grow factorially in magnitude with the cumulant order and oscillate as functions of essentially any parameter as well as of the cumulant order. We describe the theory behind this prediction which is subsequently illustrated with examples in Sec. V.

Section V is devoted to two Markovian transport models of current research interest, which we use to illustrate our recursive scheme and the generic behavior of high-order cumulants discussed in Sec. IV. We start with a model of transport through a two-level quantum dot developed by Belzig.Belzig (2005) Due to the relatively simple analytic structure of the model, it is possible to write down a closed-form expression for the CGF, allowing us to develop a thorough understanding of the behavior of high-order cumulants obtained using our recursive scheme. We study the large deviation function of the system,Touchette (2009) which describes the tails of the distribution of measurable currents, and discuss how it is related to the cumulants.

The second example concerns charge transport coupled to quantized mechanical vibrations as considered in a recent series of papers on transport through single moleculesBoese and Schoeller (2001); Braig and Flensberg (2003); McCarthy et al. (2003); Mitra et al. (2004); Koch and von Oppen (2005); Koch et al. (2005); Pistolesi et al. (2008) and other nano-electromechanical systems.Gorelik et al. (1997); Armour and MacKinnon (2002); Fedorets et al. (2002); Novotný et al. (2003, 2004); Flindt et al. (2004); Fedorets et al. (2004); Haupt et al. (2006); Rodrigues et al. (2007); Hübener and Brandes (2007); Harvey et al. (2008); Körting et al. (2009); Hübener and Brandes (2009); Cavaliere et al. (2009); Harvey et al. (2009) Due to the many oscillator states participating in transport the matrix representations of the involved operators are of large dimensions and it is necessary to resort to numerics. We demonstrate the numerical stability of our recursive algorithm up to very high cumulant orders () and show how oscillations of the cumulants can be used to extract information about the analytic structure of the cumulant generating function. We calculate the large deviation function and show that it is highly sensitive to the damping of the vibrational mode.

Section VI concerns the counting statistics of non-Markovian systems. We consider a model of non-Markovian electron transport through a Coulomb-blockade double quantum dot embedded in a dissipative heat bath and coupled to electronic leads. The dynamics of the charge populations of the double dot can be described using a non-Markovian GME whose detailed derivation is presented in App. C. We study the behavior of the first three cumulants thus extending previous studies that have been restricted to the noise.Aguado and Brandes (2004a, b) We focus in particular on the influence of decoherenceKießlich et al. (2007) on the charge transport statistics. Finally, we discuss possible subtleties associated with non-Markovian dynamics and we provide the reader with a unifying point of view on a number of results reported in previous studies as well as in the examples discussed in this paper.

Our conclusions are stated in Sec. VII. Appendix A describes the numerical algorithms used to solve the recursive equations for high-order cumulants, while Apps. B and C give detailed derivations of the Markovian GME for the vibrating molecule and the non-Markovian GME for the double-dot system, respectively.

Ii Generalized Master equation

The generic transport setup under consideration in this work is depicted in Fig. 1: A nanoscopic quantum system is connected by tunneling barriers to two electronic leads, allowing for charge and energy exchange with the leads. Typically, the quantum system consists of a discrete set of (many-body) quantum states. Moreover, the system is coupled to an external heat bath to and from which energy can flow. We will be considering a transport configuration, where a bias difference between the leads drives electrons through the system.

Figure 1: Generic transport setup. A quantum system is connected to electronic leads and a heat bath. A bias difference between the leads drives electrons through the system, which can exchange energy with the surrounding heat bath. The system is described by the -resolved density matrix (see text), where is the number of electrons that have been collected in the right lead during the time span . The probability distribution of is denoted as .

The quantum system is completely described by its (reduced) density matrix , obtained by tracing out the environmental degrees of freedoms, i. e. the electronic leads and the heat bath. It is, however, advantageous to resolve into the components , corresponding to the number of electrons that have tunneled through the system during the time span .Makhlin et al. (2001); Shelankov and Rammer (2003); Wabnig et al. (2005) The -resolved density matrix allows us to study the statistics of the number of transferred charges, similarly to well-known techniques from quantum optics.Cook (1981); Lenstra (1982); Plenio and Knight (1998) We note that the (un-resolved) density matrix can always be recovered by summing over , . For bi-directional processes, the number of tunneled electrons can be both positive and negative.

The major focus in the literature has been on systems obeying Markovian dynamics;Blum (1996); Gardiner and Zoller (2008); Alicki and Lendi (2007) however recent years have witnessed an increased interest in non-Markovian processes as well.Wilkie (2000); Aissani and Lendi (2003); Breuer et al. (2006); Bellomo et al. (2007); Breuer (2007); Budini (2008); Timm (2008); Breuer et al. (2009) In this spirit we consider a generic non-Markovian generalized master equation (GME) of the form


obtained by tracing out the electronic leads and the heat bath. An equation of this type arises for example in the partitioning scheme devised by Nakajima and Zwanzig,Zwanzig (2001) and in the real-time diagrammatic technique for the dynamics of the reduced density matrix on the Keldysh contour as described in Refs. Schoeller and Schön, 1994; König et al., 1996a, b; Makhlin et al., 2001; Braggio et al., 2006. The memory kernel accounts for the dynamics of the system taking into account the influence of the degrees of freedom that have been projected out, e. g. the electronic leads and the heat bath. Here we have assumed that the system is not explicitly driven by any time-varying fields, such that the kernel only depends on the time difference . Additionally, we assume that the number of electrons that have been collected in the right lead does not affect the system dynamics, and the kernel consequently only depends on the difference . The generic non-Markovian GME also contains the inhomogeneity which accounts for initial correlations between system and environment.Zwanzig (2001); Flindt et al. (2008) Typically, and decay on comparable time scales, and thus vanishes in the long-time limit of Eq. (1).

At this point, we note that while our method for extracting cumulants works for any GME which satisfies certain, rather general conditions specified in detail in Sec. III, the physical meaningfulness of the results nevertheless depends crucially on a consistent derivation of the -resolved memory kernel . In Sec. VI.2 we discuss various subtleties associated with a proper derivation of the memory kernel for non-Markovian systems. In this work we use the notion of the Markovian limit of a general non-Markovian GME in a somewhat loose manner, namely by referring to the “Markovian” limit of Eq. (1) as the case, where and . We use this terminology for the ease of notation, although we are aware that the proper Markovian limit under certain circumstances may actually be different. For an example of this, we refer the reader to Ref. Dümcke and Spohn, 1979, where it is demonstrated that the correct Markovian limit for weak coupling theories should be performed in the interaction picture. Since this procedure only influences the off-diagonal elements in the weak coupling regime, we ignore this subtlety in the rest of the paper as we will not be considering such cases. In relevant situations this difference should be taken into account — it would, however, only lead to a reinterpretation of the non-Markovian corrections studied in Sec. VI.2. The issue of non-Markovian behavior, its nature and distinction from Markovian approximations, is a nontrivial and timely topicWolf et al. (2008); Breuer et al. (2009); Rivas et al. (2009) which we only touch upon briefly in this work, but our formalism paves the way for systematic studies of such problems in the context of electronic noise and counting statistics, for example as in Ref. Zedler et al., 2009.

ii.1 Counting statistics

In the following we introduce the notion of cumulants of the charge transfer probability distribution, and derive a formal expression for the cumulant generating function CGF from the GME (1). The probability distribution for the number of transferred particles is obtained from the -resolved density by tracing over the system degrees of freedom,


Obviously, probability must be conserved, such that . In order to study the cumulants of it is convenient to introduce a cumulant generating function (CGF) via the definition


from which the cumulants follow as derivatives with respect to the counting field at ,


Alternatively, one can write


which defines the -dependent density matrix


By going to Laplace space via the transformation


Equation (1) transforms to an algebraic equation reading


This equation can be solved formally by introducing the resolvent


and writing


Finally, inverting the Laplace transform using the Bromwich integral we obtain for the CGFFlindt et al. (2008)


where is larger than the real parts of all singularities of the integrand.

Equation (11) is a powerful formal result for the CGF, and, as we shall see, it also leads to practical schemes for calculating current fluctuations. In this work, we concentrate on the zero-frequency cumulants, determined by the long-time limit of the CGF. The case of finite-frequency noise,Flindt et al. (2008) where the inhomogeneity plays an important role, will be considered in future works.

ii.2 Notational details

Throughout this paper we will use the superoperator notation previously described in Ref. Flindt et al., 2004 and also used in a number of other works.Flindt et al. (2005a, b); Hübener and Brandes (2007); Flindt et al. (2008); Brandes (2008); Harvey et al. (2008); Körting et al. (2009); Hübener and Brandes (2009); Emary (2009); Emary et al. (2009); Schaller et al. (2009); Harvey et al. (2009); Urban and König (2009); Lindebaum et al. (2009); Zhong and Cao (2009); Wu and Timm (2010); Domínguez et al. (2010) Using this notation, standard linear algebra operations can conveniently be performed, analytically and numerically. Within the formalism, the memory kernel , the resolvent , and other operators that act linearly on density matrices, are referred to as superoperators and denoted by calligraphic characters. Conventional quantum mechanical operators, like the density matrix , acting in the conventional quantum mechanical Hilbert space, can be considered themselves to span a Hilbert space, referred to as the superspace. The superoperators act in the superspace, while conventional quantum mechanical operators are considered as vectors using a bra(c)ket notation, i.e., , where is a conventional quantum mechanical operator, and is the corresponding ket in the superspace. Double angle brackets are used here in order to avoid confusion with conventional kets. In numerical calculations, bras and kets are represented by vectors, while superoperators are represented by matrices. The inner product between bras and kets is defined as . Since the involved superoperators, like and , are not hermitian, their eigenvalues are generally complex. In such cases, left and right eigenvectors corresponding to a particular eigenvalue are not related by hermitian conjugation. The left eigenvector, or bra, corresponding to an eigenvalue is therefore denoted with a tilde, e.g. , to avoid confusion with the hermitian conjugate of the corresponding right eigenvector, or ket, .

Iii Zero-frequency current cumulants

In this section we derive the recursive method for evaluating the zero-frequency current cumulants. We first define the zero-frequency cumulants of the current as


where . As we shall show below, the cumulants of the passed charge become linear in at long times such that , and the zero-frequency current cumulants are thus intensive quantities (with respect to time). Thus, in the long-time limit , and we use these two normalized quantities interchangeably throughout the paper.

In order to find the long-time limit of the CGF, we consider the formal solution Eq. (11). The memory kernel is assumed to have a single isolated eigenvalue , which for is zero, corresponding to the stationary limit of , i.e., for large . Here, is the normalized solution to . We exclude cases, where the zero-eigenvalue is degenerate due to two or more uncoupled sub-systems.van Kampen (2007) In the bracket notation is denoted as . The corresponding left eigenvector can be found be noting that the memory kernel with conserves probability for any . This can be inferred from the GME in Laplace space: For normalized density matrices with , we have , and Eq. (8) yields


It is generally possible to choose an initial state such that . The kernel does not depend on the choice of initial state and since Eq. (13) holds for any normalized density matrix we deduce that . In the bracket notation this equality can be expressed as with the left eigenvector in the superspace corresponding to the identity operator in the conventional Hilbert space. This moreover implies thatBraggio et al. (2006)


We next examine the eigenvalue which we assume evolves adiabatically from with small and . It is convenient to introduce the mutually orthogonal projectors




with developing adiabatically from for small and . Here, and are the left and right eigenvectors corresponding to , which develop adiabatically from and , respectively. In terms of and the memory kernel can be partitioned as


In deriving this expression we used


Using the partitioning, Eq. (17), the resolvent becomes


For the first term of the resolvent has a simple pole at , which determines the long-time limit, i.e. it corresponds to the stationary state . We denote the pole at by . All singularities of the second term have negative real parts and do not contribute in the long-time limit. Again, we assume adiabatic evolution of the pole from with small , such that is the particular pole that solvesBraggio et al. (2006); Flindt et al. (2008)


With small , the other singularities still have more negative real parts and the pole again determines the long-time behavior. From Eq. (11) we then find for large




From the definition of the zero-frequency current cumulants in Eq. (12) we then establish that


We note that the CGF in the long-time limit and thus the zero-frequency cumulants do not depend on the initial state and the inhomogeneity . In contrast, both and must be appropriately incorporated in order to calculate the finite-frequency noise.Flindt et al. (2008) Equations (20) and (23) form the main theoretical result of this section, generalizing earlier results for Markovian systems.Bagrets and Nazarov (2003); Flindt et al. (2005a) In the Markovian limit, the memory kernel and the corresponding eigenvalue close to 0 have no -dependence, and Eq. (20) immediately yields ,Bagrets and Nazarov (2003); Flindt et al. (2005a) where is the eigenvalue of the -independent kernel, which goes to zero with going to zero, i.e. .

Although, we have formally derived an expression for the CGF, it may in practice, given a specific memory kernel , be difficult to determine the eigenvalue including its dependence on and . Moreover, the solution of Eq. (20) itself poses an additional problem, which needs to be addressed in the non-Markovian case. In the Markovian limit, only derivatives of the eigenvalue with respect to the counting field need to be determined. However, with the superoperator being represented by a matrix of size , there is no closed-form expression for the eigenvalue already with . The immediate alternative strategy would then be to calculate numerically the eigenvalue and the derivatives with respect to . Typically, however, this is a numerically unstable procedure, which is limited to the first few derivatives.Press et al. (2007) Consequently, we devote the rest of this section to the development of a numerically stable, recursive scheme that solves Eqs. (20) and (23) for high orders of cumulants, including in the non-Markovian case.

iii.1 Recursive scheme

iii.1.1 The Markovian case

We consider first the Markovian case,Flindt et al. (2005a); Baiesi et al. (2009) before proceeding with the general non-Markovian case. In the Markovian case, the memory kernel has no -dependence, and the current cumulants are determined by the eigenvalue which solves the eigenvalue problem


where . We find the eigenvalue using perturbation theory in the counting field in a spirit similar to that of standard Rayleigh-Schrödinger perturbation theory. To this end we introduce the unperturbed operator and the perturbation such that


We can then write


where we have used and chosen the conventional normalization . We moreover employ the shorthand notation and for the projectors introduced in the previous section, and write


consistently with the choice of normalization. Using that , Eq. (24) can be written


Next, we introduce the pseudo-inverseFlindt et al. (2004, 2005a) defined as111We note that our definition of the pseudo-inverse differs from the Moore-Penrose pseudo-inverse , which is implemented in many numerical software packages, e.g., in Matlab. However, when projected on the regular subspace by it reduces to our pseudo-inverse, i.e. .


The pseudo-inverse is a well-defined object, since the inversion is performed in the subspace corresponding to , where is regular. By applying on both sides of Eq. (28) we find


which combined with Eq. (27) yields


Equations (26) and (31) form the basis of the recursive scheme developed below.

We first Taylor expand the eigenvalue , the eigenvector , and the perturbation , around as


where we have used that and . Inserting these expansions into Eqs. (26) and (31), and collecting terms to same order in , we arrive at a recursive scheme reading


for . The recursive scheme allows for systematic calculations of cumulants of high orders.

As illustrative examples we evaluate the first three current cumulants using the recursive scheme,


having used and , since . The subscript reminds us that these results hold for the Markovian case. The expressions (34) for the first three cumulants are equivalent to the ones derived in Ref. Flindt et al., 2005a, albeit using a slightly different notation. Importantly, the recursive scheme presented here allows for an easy generation of higher order cumulants, either analytically or numerically.

iii.1.2 The non-Markovian case

We now proceed with the non-Markovian case, where we first need to consider the eigenvalue problem


where is the particular eigenvalue for which . The basic equations, Eqs. (26) and (31), are still valid, provided that , , , and are replaced by , , , and , respectively, i.e.,




Again, we Taylor expand all objects around , but in this case also around ,


with by definition and , since . Inserting these expansions into Eqs. (36) and (37) and collecting terms to same orders in and , we find a recursive scheme reading


In case the memory kernel has no -dependence, corresponding to the Markovian case, only terms with are non-zero, and the recursive scheme reduces to the one given in Eq. (33). In particular, the coefficients equal the current cumulants in the Markovian limit of the kernel, .

In the non-Markovian case, we need to proceed with the solution of Eq. (20) for and extract the current cumulants . Inserting the expression for in Eq. (23) into Eq. (20) and using the expansion of given in Eq. (38), we find


Collecting terms to same order in , we find


in terms of the auxiliary quantity


where only terms in the sums for which should be included. For , we have . The auxiliary quantity can also be evaluated recursively by noting that


with the boundary conditions , , and .

When combined, Eqs. (39, 41, 43) constitute a recursive scheme which allows for numerical or analytic calculations of cumulants of high orders in the general non-Markovian case. As simple examples, we show the first three cumulantsFlindt et al. (2008) obtained from Eqs. (41), (43), in terms of the coefficients


In general, the ’th current cumulant contains the coefficients


with . However, coefficients of the form are zero since as discussed below Eq. (13) and it thus suffices to consider . From Eq. (39) it follows that depend only on with and so that we can conclude that the ’th cumulant of the current depends at maximum on the ’th time-moment of the memory kernel . In particular, this implies that the mean current is a purely Markovian quantity depending only on the time-integrated memory kernel while the second and higher order cumulants deviate from the results in the Markovian case.Braggio et al. (2006)

The coefficients can be found from Eq. (39). Coefficients of the form only contain zeroth order terms in and are, as already mentioned, equal to the current cumulants in the Markovian limit, i.e.,


For the other coefficients entering the expressions in Eq. (44) for the first three non-Markovian current cumulants, we find


Again, as in the Markovian case, higher order cumulants including the coefficients are readily generated, analytically or numerically. The results presented here can be generalized to the statistics of several different counted quantities as in Ref. Sánchez et al., 2007; Sánchez et al., 2008a, and cross-correlations can be evaluated using the same compact notation developed in this work.Braggio et al. (2009a)

iii.2 Notes on evaluation

As previously mentioned, the size of the memory kernel could in practice hinder the calculation of and the solution of Eq. (20), and thus the evaluation of the current cumulants. The recursive scheme described above, however, only relies on the ability to solve matrix equations and perform matrix multiplications. Both of these operations are numerically feasible and stable, even when the involved matrices are of large dimensions. In general, the recursive scheme requires the following steps: The stationary state must be found by solving


with the normalization requirement . Secondly, the and derivatives of the memory kernel must be found


for . Typically, the dependence on the counting field enters matrix elements in an exponential function (see e. g. Refs. Nazarov, 2003, Bagrets and Nazarov, 2003, Braggio et al., 2006 and examples in Secs. V and VI), e. g. as a factor of , for which the derivatives with respect to are easily found analytically. The -dependence of the matrix element can be written


such that


The integration over time can be performed in a numerically stable manner for arbitrary ,Press et al. (2007) thereby avoiding taking numerical derivatives with respect to .

Finally, matrix multiplications have to be performed. Here, special attention has to be paid to terms involving the pseudo-inverse , i.e. , where for example has the form in the expression for the coefficient in Eq. (47). In order to evaluate such expressions we introduce as the solution(s) to


such that


which can be verified by applying on both sides of Eq. (52) and using that and . The projector in Eq. (52) ensures that the right hand side lies in the range of , and since is singular, the equation has infinitely many solutions. The solutions can be written


where is a particular solution to Eq. (52), which can be found numerically. We then obtain by applying to according to Eq. (53) and find


since .

In App. A we describe a simple numerical algorithm for solving Eqs. (48) and (52). For very large dimensions of the involved matrices, it may be necessary to invoke more advanced numerical methods to solve these equations.Flindt et al. (2004) Numerically, the recursive scheme is stable for very high orders of cumulants (up to order ), which we have tested on simple models. The results presented in this work have all been obtained using standard numerical methods as the one described in App. A.

Iv Asymptotics of high-order cumulants

Before illustrating our methods in terms of specific examples, we discuss the asymptotic behavior of high-order cumulants. As some of us have recently shown certain ubiquitous features are expected for the high-order cumulants.Flindt et al. (2009) In particular, the absolute values of the high-order cumulants are expected to grow factorially with the cumulant order. Moreover, the high-order cumulants are predicted to oscillate as functions of basically any parameter, as well as of the cumulant order. This behavior was confirmed experimentally by measurements of the high-order transient cumulants of electron transport through a quantum dot.Flindt et al. (2009) In the experiment, the transient cumulants indeed grew factorially with the cumulant order and oscillated as functions of time (before reaching the long-time limit), in agreement with the general prediction. For completeness, we repeat here the essentials of the theory underlying these asymptotic properties of high-order cumulants.

The asymptotic behavior of high-order cumulants follows from straightforward considerations. In the following we denote the CGF by , where represents the set of all parameters needed to specify the system; whether the dynamics is Markovian or non-Markovian is irrelevant. In general, we can assume that the CGF has a number of singularities in the complex- plane at , , which can be either poles or branch-points. Typically, the positions of the singularities depend on . Exceptions, where the CGF has no singularities, do exist, e. g. the Poisson process, whose CGF is given by an exponential function, but we exclude such cases in the following.

Close to a singularity , we can write the CGF as


for some and , determined by the nature of the singularity. For example, for a finite-order pole denotes the order of the pole, while would correspond to the branch point of a square-root function. Logarithmic singularities can be treated on a similar footing with only slight modifications.Flindt et al. (2009) The derivatives with respect to the counting field are now




for . As the order is increased this approximation becomes better away from the singularity at according to the Darboux theorem.Dingle (1973); Berry (2005); Flindt et al. (2009) For sufficiently high , the cumulants of the passed charge can thus be written


where the sum runs over all singularities of the CGF. Here, we have written the singularities as


where is the modulus of the singularity and is the corresponding complex argument. In general, the singularities together with the factors come in complex conjugate pairs, ensuring that the expression in Eq. (59) is real.

From Eq. (59) we deduce that the cumulants grow factorially in magnitude with the order due to the factors given in Eq. (58). We also see that the high-order cumulants are determined primarily by the singularities closest to zero. Contributions from other singularities are suppressed with the relative distance from zero and the order , and can thus be neglected for large . Importantly, we observe that the high-order cumulants become oscillatory functions of any parameter among that changes as well as of the cumulant order [see also Eq. (61) below]. We refer to these ubiquitous features, which should occur in a large class of transport processes, as universal oscillations. For example, we expect oscillations of high-order cumulants for basically any transport process described by a GME, since the CGF for these systems typically have logarithmic singularities at finite timesFlindt et al. (2009) or square-root branch points in the long-time limit.Bender and Orszag (1999) Factorial growth and oscillations as functions of various parameters can be found in several independent studies of high-order cumulants,Pilgram and Büttiker (2003); Förster et al. (2005, 2007); Flindt et al. (2008); Urban et al. (2008); Khoruzhenko et al. (2009); Prolhac and Mallick (2009); Golubev et al. (2010); Hassler et al. (2010) as well as in the recent experiment described in Ref. Flindt et al., 2009, demonstrating the generality of the phenomenon. Similar observations and discussions can also be found in quantum opticsDodonov et al. (1994) and high-energy physics,Dremin and Hwa (1994); Bhalerao et al. (2003, 2004) further confirming the prediction. We note that in the long-time limit, the positions of the dominating singularities are no longer time dependent,Flindt et al. (2009) and the cumulants cease to oscillate as functions of time. Instead, the cumulants of the passed charge become linear in time, as previously discussed in Sec. III.

A simple (and common) situation arises if only two complex conjugate singularities, and , are closest to zero. In that case, Equation (59) immediately yields


Using this expression we can determine the positions of the dominating singularities from numerical calculations of the high-order cumulants as we shall demonstrate in the second example considered in Sec. V. We note that while the factorial growth and the oscillations are system independent, other features, for example the frequency of the oscillations, are determined by the particular details of the system under consideration.

Finally, we mention the Perron-Frobenius theorem regarding stochastic matricesBaiesi et al. (2009); Demboo and Zeitouni (1998) which implies that the CGFs considered in Sec. V must be analytical functions at least in a strip along the real axis in the complex- plane. This has important consequences especially for the nature of the high-order cumulants which rests heavily on the analytical properties of the CGF. We illustrate this statement in both examples in Sec. V.

V Markovian Systems

v.1 Electron bunching in a two-level quantum dot

In our first example we study electron bunching in transport through a two-level quantum dot as described by Belzig in Ref. Belzig, 2005. Due to the relatively simple analytical structure of the model, it is possible to illustrate the concepts of universal oscillations introduced above. The model allows us to test the accuracy of our numerical calculations of high-order cumulants against analytic expressions.

We start by summarizing the setup in Belzig’s model. Consider a single quantum dot with two single-particle levels coupled to voltage-biased source and drain electrodes. The two levels serve as parallel transport channels. Due to strong Coulomb interactions on the quantum dot only one of the levels can be occupied at a time. The system exhibits super-Poissonian bunching transport in cases where both levels are coupled by the same rate to the, say, left lead, whose Fermi level is kept well above both levels, while the couplings to the other lead are markedly different, such that one level is coupled to the right lead by the rate and the other by with . This situation can arise for example, if the two levels are situated above and slightly below, respectively, the Fermi level of the right lead at a finite electron temperature.

This particular configuration leads to bunching of electrons in the transport due to the existence of the blocking state: if the dot is empty there is equal probability for either of the two levels to be filled. Current runs easily through the first level, while the other level effectively is blocked, or more precisely, the transport through the level is limited by the very small right rate , constituting a bottleneck. The transport thus proceeds in bunches of electrons passing intermittently through the first level separated by quiet periods of blocked transport when the other level is occupied. This bunching effects leads to super-Poissonian noise with a Fano factor above unity. For more detailed discussions of the model as well as its generalizations to many levels, the reader is referred to Ref. Belzig, 2005.

The counting statistics of the system can be obtained from a Markovian rate equation for the probability vector , containing the (-resolved) probabilities for the quantum dot to be empty, or the first (, non-blocking) or second (, blocking) level being occupied, respectively. The corresponding -dependent rate matrix reads


Here, we have rescaled the time and set while renaming in order to simplify the analytic results in the following. We have also made a minor modification of the model in Ref. Belzig, 2005 by including the back-flow into the blocking level from the right lead. This modification, however, changes only slightly the detailed quantitative results, while leaving the main qualitative features identical in the limit of interest .

Since the model involves only three states, the CGF can be found analytically in the long-time limit. The full expression is too lengthy to be presented here, but in the limit , it reduces to the result by BelzigBelzig (2005) (also for our slightly modified model; note, however, the opposite sign convention for the CGF in Ref. Belzig, 2005)


Clearly, the CGF has simple poles at


with the pole being closest to 0. However, according to the Perron-Frobenius theorem mentioned in Sec. IV the CGF cannot have singularities on the real -axis.

In order to illustrate this point, we consider the expected behavior of the high-order cumulants based on the CGF above. Close to the singularity , we approximate the CGF by the first non-zero term of the Laurent series


This corresponds to Eq. (56) with and . From Eq. (59) we then obtain a simple asymptotic expression for the high-order cumulants reading


Here, the subscript indicates that the expression has been obtained using the approximate CGF in Eq. (63) with only a single singularity closest to zero. In Table 1 we compare the asymptotic expression with results for the first six cumulants obtained by direct differentiation of the CGF in Eq. (63). The asymptotic results are very close to the exact derivatives of the approximate CGF. Despite the good agreement with the approximate results, the asymptotic expression in Eq. (66) does not reproduce our numerically exact results, also shown in the table, obtained using our recursive scheme. In particular for high orders, the asymptotic expression starts to deviate significantly from the numerically exact results.

2 3 4 5 6
Single-pole approx. 2.000 6.000 26.00 150.0 1082 9366
Single-pole asympt. 2.081 6.006 25.99 150.0 1082 9366
Numerics 1.978 5.880 25.18 143.3 1017 8644
Table 1: Normalized zero-frequency current cumulants for transport through a two-level quantum dot. Single-pole approximation results have been obtained by direct differentiation of the CGF in Eq. (63) or its asymptotic expression Eq. (66), respectively. The numerically exact results have been obtained using our recursive scheme and the rate matrix in Eq. (62) with and .
Figure 2: High-order cumulants and large deviation function for bunching transport through a two-level quantum dot. Left and central panels show comparisons between exact numerics and the single pole approximation stemming from Eq. (63) for two different values of and . The asymptotic expression in Eq. (61) based on a pair of complex conjugate singularities is shown with full lines. Notice that . The right panel shows a comparison of the large deviation function (LDF) obtained from exact numerics and the single pole approximation in Eq. (71), respectively.

As anticipated above, these deviations can be traced back to the expression in Eq. (63), that we obtained in the limit . In order to proceed from here, we return to the full expression for the CGF in the long-time limit (not shown). A careful analysis reveals that, in fact, there is a pair of complex conjugate singularities closest to zero, and not just a single pole. The two singularities, denoted as and , correspond to branch points of a square-root, and for small the position of the branch point is


Clearly, for small the branch points are close to the position of the single pole . However, for any finite , the two branch points have small, but finite, imaginary parts thus complying with the Perron-Frobenius theorem. The singularity structure around the branch point is characterized by Eq. (56) with and , and we can then use the asymptotic expressions in Eq. (61) for the high-order cumulants. In the left and central panels of Fig. 2 we compare this expression, and the single-pole approximation in Eq. (66), with numerically exact results obtained using our recursive scheme for and two different values of .

Figure 2 shows several important features. Firstly, the (scaled) high-order cumulants indeed behave in an oscillatory manner as function of the cumulant order , which coincides with the cosine part of Eq. (61). Obviously, for smaller the period of the oscillations, determined by , is longer in accordance with Eq. (61). Furthermore, for the small value of , the asymptotic form of the high-order cumulants is reached around , while the single-pole approximation agrees well for lower orders, . For the higher value of , significant deviations from the single-pole behavior begin already for the fourth cumulant, while the asymptotic oscillatory form holds from around . Notice the importance of the exact analytical knowledge of the singularities — even though (together with ) may seem a very small number justifying the usage of the single pole approximation, we see from Eq. (67) that the imaginary part of the pole and its argument scale like , thus invalidating the single-pole approximation far earlier than expected from a linear-in- scaling assumption.

A complementary view on the charge transport statistics is provided by the large deviation function (LDF),Touchette (2009) which quantifies deviations of measurable currents from the average value. The LDF is obtained from the probability distribution


and is defined as the long-time limit of , where is the current. For long times, we have and the integral can be evaluated in the saddle-point approximation with the saddle-point given by the solution to the saddle-point equation


The saddle-point equation implies a parametric dependence of the saddle-point on the current . Using the saddle-point approximation, the LDF becomes


We first solve the saddle-point equation for the approximate CGF in Eq. (63) and find


where and the subscript again reminds us that the expression has been obtained using the approximate CGF with only a single singularity closest to zero. Obviously, the current must be positive (), since transport is unidirectional.

Also for the LDF, we can compare the analytic approximation with numerical exact results. To this end, we need to solve the saddle-point equation taking as starting point the kernel in Eq. (62). The derivative of the eigenvalue is now calculated using the Hellman-Feynman theorem, writing