Finite-frequency counting statistics of electron transport: Markovian Theory

We present a theory of frequency-dependent counting statistics of electron transport through nanostructures within the framework of Markovian quantum master equations. Our method allows the calculation of finite-frequency current cumulants of arbitrary order, as we explicitly show for the second- and third-order cumulants. Our formulae generalize previous zero-frequency expressions in the literature and can be viewed as an extension of MacDonald's formula beyond shot noise. When combined with an appropriate treatment of tunneling, using, e.g. Liouvillian perturbation theory in Laplace space, our method can deal with arbitrary bias voltages and frequencies, as we illustrate with the paradigmatic example of transport through a single resonant level model. We discuss various interesting limits, including the recovery of the fluctuation-dissipation theorem near linear response, as well as some drawbacks inherent of the Markovian description arising from the neglect of quantum fluctuations.


Introduction.
Transport of electrons through nanoscopic conductors is a very powerful tool to learn about interactions and to characterize quantum systems [1]. Examples include the quantum Hall effect [2], weak localization [3], or universal conductance fluctuations [4]. Transport processes are governed by tunneling events, which are stochastic in nature. It is therefore natural to expect that the statistics of these tunneling events will be strongly influenced by interactions and quantum effects. Interestingly, these statistics, which can be analyzed by studying current fluctuations, contain a great deal of new information beyond that provided by dc transport [5,6]. In particular, the second-order current correlation function (noise), can be used to determine the effective charge [7] and the statistics of the quasiparticles [8], and to reveal information on the transmission properties of the conductor [5]. Moreover, current correlations can be used to learn about entanglement [9,10], quantum coherence [11], and the deep connection that exists with fluctuation theorems [12,13,14]. Further information can be gained from noise at high frequencies which is valuable for extracting internal energy scales in systems such as quantum dots [15], spin-valves [16], Cooper Pair Boxes [17], diffusive wires [18] or chaotic cavities [19].
A proper treatment of fluctuations in non-equilibrium transport is needed to address the problems listed above. While the list of theoretical available is too large to be given here, it is safe to say that they can be divided roughly into three families: the scattering approach [39,5], the Keldysh Green's functions method [40], and the various quantum master equation (QME) treatments [41] (for a recent overview of transport in this last context see e.g. Ref [42]). A theory of counting statistics of electron transport was first formulated by Levitov and Lesovik for noninteracting electrons using the scattering formalism [43], and later works enabled the treatment interacting problems [44,45]. QME approaches followed [46], and prooved particularly useful for studying systems in the Coulomb Blockade regime. Recent avances within this last scheme also involve studies of the counting statistics including non-Markovian dynamics [47,48,49,50].
In Ref. [51] we developed a method for calculating high-order current correlations at finite frequencies in the context of Markovian QMEs, and the aim of this paper is to significantly extend that work. We provide a detailed derivation of our multitime generating function, Eq. (27). We present a new approach to derive finitefrequency cumulants from this expression, noticing that only the stationary solution of the problem is required. We give analytical expressions for the second and thirdorder current cumulants (noise and skewness respectively). These results generalize previous zero-frequency expressions in the literature and recover the finite frequency shot noise expressions [52,53] obtained using the MacDonald's formula [54]. Our method can thus be viewed as a generalization of this formula, as it allows us to obtain high-order spectra such as the frequency-dependent skewness (Eq. (35)). To illustrate the formalism we study the case of a single resonant level (SRL) model, and compare it with the exact solution and with the non-equilibrium version of the fluctuation-dissipation theorem, derived in various approaches, such as for tunnel junctions, or for the weak cotunneling regime in quantum dots [55,56,57,58].
The paper is organized as follows. In section 2, we present our formalism of finite-frequency cumulants in the context of Markovian QMEs. Subsections 2.1 and 2.2 are devoted to establish the general framework of full counting statistics. In subsection 2.3 we derive a multitime cumulant generating function. Subsection 2.4 shows how to obtain finite-frequency cumulants of the current distribution. Here, exact equations for the frequency-dependent noise and skewness are given. We end this part with special emphasis on how to calculate the counting statistics of the "total" and "accumulated" currents (subsection 2.5). We explicitly show how both current correlations and charge correlations can be calculated. In section 3, we study the example of a SRL model, providing spectra for the frequency-dependent noise and skewness and a detailed comparison with the fluctuation-dissipation theorem, the finite-frequency version of the non-equilibrium fluctuation-dissipation theorem, and the exact solution of the SRL model. First we focus on the zero-frequency case (section 3.1), where the general behaviour of noise and skewness is presented as a function of different system parameters. In section 3.2 we extend this study to the finite-frequency case. Interestingly, we show that, even though the theory does not contain quantum fluctuations, the Markovian limit is basically exact in transport configurations, level within the bias voltage window, as long as ω ≫ eV or ω ≪ eV . In intermediate situations, where ω ∼ eV , or with the level outside the bias windows, the Markovian limit fails at finite frequencies due to its lack of quantum fluctuations. We also demonstrate that the noise spectra for particle currents and the ones for total currents significantly deviate from each other, even for large asymmetric coupling to the leads, namely Γ R /Γ L = 1. In section 4 we summarize our results. Most of the technical details and intermediate steps of the derivations in section 2 are discussed in detail in Appendix A, where we also present a diagrammatic technique to arrive to the expressions for the cumulants shown in sec. 2.4. In section Appendix B we describe how to calculate the kernel of the QME to lowest (sequential) order using perturbation theory in the Liouville space.

Quantum master equation
We are interested in phenomena which can be described by an equation of the typė where L is the so called Liouvillian, governing the evolution, and ρ is the density operator of the total system. Specifically, our theory will be useful to processes amenable to the counting of a classical stochastic variable n, which can be, for example, the number of particles that have undergone a particular process in the system. We will focus in particular in transport systems, consisting on a central region, with a known set of many-body eigenstates {|a }, and respective eigenenergies {E a }, attached to non-interacting electronic leads at different chemical potentials. This set-up can be described by a Hamiltonian that takes the form H = H S + H R + H T , where H S and H R refer to system and leads respectively, and H T is the coupling between them. The different terms can be written as Here c † ηα creates an electron with quantum numbers η in lead α, and d m annihilates an electron from site m in the central region. ε ηα are the eigenenergies of the electrons in the lead α, V ηαm is the coupling energy between a state in contact α and the level m in the system. µ α the chemical potential of lead α, and that allows the system to be driven out of equilibrium.
Given the Hamiltonian, the full evolution can be written in terms of the Liouvillian by using the von Neumann equatioṅ We are actually interested in the dynamics concerning the central system. We therefore trace out the reservoir degrees of freedom, arriving at an equation of motion for the reduced density operator ρ S (t) ≡ Tr R {ρ(t)}, that, if we assume the dynamics of the reservoirs to be much faster than that of the central system, can be approximated as a Markovian master equatioṅ where W is a kernel driving the dynamics of the system. The charge flow through the conductor is governed by the stochastic hopping of electrons in and out of the central region. This processes are susceptible to classical counting and thus the reduced density operator can be unravelled into components ρ S (n α , t), corresponding to having n α = ..., −1, 0, 1, ... extra electrons in lead α [59,60]. The kernel can also be split as where W α ± refers to the physical process in which one electron in created (+) or annihilated (-) at lead α, and W 0 corresponds to the part in which no tunneling processes take place in that lead. It can be shown that ρ S (n α , t) fulfills the equation (see e.g. [61]): valid provided that only single particle tunneling processes occur. Although Eq. (7) focuses on the single counting at a particular lead α, it can be generalized to account for tunneling processes of k particles at the different system-reservoir junctions: ρ S (n 1 , . . . , n M , t) = W 0 ρ S (n 1 , . . . , n M , t) where W 0 is the part in which the number of particles is not changed in the central region, M the number of leads, and k labels the process in which k particles "jump" at a time. Unfortunately, solving equation (8) in the n-space requires truncation to a certain n and diagonalization of a tridiagonal matrix. It is therefore more convenient to solve it taking its Fourier transform. Multiplying (8) by e in1χ1 . . . e inM χM and summing over n 1 , . . . , n M we obtaiṅ Here the counting field χ refers implicitly to all counting fields, and we have defined ρ S (χ, t) ≡ n1,...,nM e in1χ1 . . . e inM χM ρ S (n 1 , . . . , n M , t) and W(χ) ≡ W 0 + α,± W α ± e ±iχα . For a time-independent kernel, the solution to Eq. (9) is with the time evolution operator Ω(χ, t − t 0 ) := e W(χ)(t−t0) .

Full Counting Statistics
Importantly, the knowledge of the system's density operator resolved in n allows us to obtain the the full counting statistics (FCS) of the system, that is the probability distribution P (n 1 , . . . , n M , t) of the number of electrons transmitted through the system-lead junctions. This is accomplished by noting that P (n 1 , . . . , n M , t) = Tr{ρ S (n 1 , . . . , n M , t)}.
Transforming the probability distribution to the χ-space, we have the moment generating function (MGF) where now n refers implicitly to n 1 , . . . , n M . Using Eq. (10) one has equation that was already derived by Bagrets and Nazarov [46]. The N -th. derivative of the MGF with respect to χ gives the N -th moment of the distribution of the number of particles that have gone in or out a particular lead α: When equation (13) is used, averages with respect to the stationary state are established by taking ρ S (t 0 ) = ρ stat S (defined by Wρ stat S = 0). This means that counting will start at a time t 0 in which the system has reached its steady state, and therefore the fluctuations we study are around this state †.
The moments of the current distribution can be calculated as ‡ This relation is important as it relates the stochastic variable n with the current of particles flowing through the system. Even though the current studied here is a classical variable, it contains quantum effects present in the system. In the formalism, these are inherited from the Liouvillian operator in equation (1). Generally, we are interested in the cumulants, rather than the moments, of the current distribution. These can be obtained from the derivatives of the cumulant generating function (CGF), defined by F (χ, t) := lnG(χ, t). Therefore we have † In the following, all the averages will be taken with respect to this steady state. An average in the Liouville space will be therefore written as is the normalized stationary system density matrix (written as a vector), and t T ≡ 0 | is the transposed trace vector that sums over the population degrees of freedom. ‡ Throughout the paper we will use e (electron charge) = k (Boltzman's constant) = (Planck's constant/2π) = 1.
where . . . c denotes cumulant average [62] and the limit t → ∞ ensures that the average is performed in the stationary state. Also, notice that the probability distribution itself can be obtained by inverse Fourier transform of the MGF.
From the χ-independent kernel of the reduced QME (6), the χ-dependence leading to equation (9) can be actually introduced in a simpler way than resolving the density operator in n and taking the Fourier transform. As we describe in Appendix B, it is enough to include counting fields in the appropriate tunneling terms of the Kernel [63], and this procedure is fully equivalent to solving a generalized Von Neumann equation: in which the time evolution in the forward (+) and backward (-) Keldysh part of the real time axis is governed by different Hamiltonians [55], specifically H ± T = η,α,m V ηαm e ±iχ/2 c † ηα d m + H.c.. Tracing out the reservoir degrees of freedom in equation (17), one can get equation (9) and proceed to obtain the FCS of the system.

Finite-frequency full counting statistics
In this paper we want to study correlations at finite frequencies, for which the scheme presented above has to be generalized. To this end, one needs to consider a joint probability distribution, P (n 1 , t 1 ; . . . ; n N , t N ) defined as the probability that n 1 electrons have undergone a particular process after a time t 1 , n 2 electrons after a time t 2 , etc. Here we focus for simplicity on a particular lead, and denote n as the number of particles transferred to (from) it, with associated counting field χ (−χ). Being straightforward to include processes at different leads.
The connection between this joint probability and the density operator (analog to Eq. (11)) is not straightforward. To connect them we first need to specify a prescription for the symmetrization of the cumulants and the probability distribution. This prescription actually depends on the detection scheme [64,65,66]. Here we assume that "classical" detection is carried out, so the detector is incapable of distinguishing emission from absorption. This means that the results we will present correspond to the fully-symmetrized version of the power spectrum where T S is the symmetrization operator, that sums over all possible time (or frequency) switchings, that is, we have for example T S I(t 1 )I(t 2 ) = I(t 1 )I(t 2 ) + I(t 2 )I(t 1 ) . The spectrum (18) can be derived from a N -time (symmetrized) CGF F (N ) , defined by where χ := (χ 1 , . . . , χ N ) T , t := (t 1 , . . . , t N ) T , n := (n 1 , . . . , n N ) T , P (N ) refers to the symmetrized joint probability, G (N ) to the multitime MGF, and the subscript T to the transpose of a column vector. That is, we have And using the property of the Fourier transform of a derivative we get Both the probability P (N ) [n, t] and the CGF F (N ) [χ, t] can be calculated from the density operator and the Kernel W if we use the Markovian approximation in the coupling with central system-reservoir coupling. Within that limit we have the evolution local in time given by (10) and also the factorization property where the symbol > constraints the times to t k > t k−1 . Notice that as we are considering the totally symmetric correlation function, we need to consider where T is the time-ordering operator [51]. P (n, t|n ′ , t ′ ) is the conditional probability of counting n electrons after time t provided that we counted n ′ electrons after time t ′ , and can be computed as where the normalization in the denominator accounts for the collapse of the state due to the measurement, as given by the von Neumann's projection postulate [67]. Ω(n, t) is the propagator in the n-space, that is, and can be extracted from equation (8) or by inverse Fourier transform of the propagator in the χ-space: An expression for the joint probability distribution in terms of propagators can be then derived using (21) together with (22) and (23). Alternatively, it can be obtained using the Chapman-Kolmogorov property for Markovian evolutions, from which we have where ν k := n k+1 −n k , τ k := t k+1 −t k and • := Tr{•ρ stat S }. Transforming expression (26) to the χ-space, we find the CGF (27) is encountered in many branches of physics such as statistical physics and field theory, where connected correlation functions are obtained by taking derivatives of the logarithm of the corresponding generating functional (the partition function, the S-matrix, etc). Note in particular the analogy with the partition function presented in Ref. [62].

Finite-frequency cumulants
Equation (27) permits us to obtain frequency-dependent current cumulants to arbitrary order. This was used in Ref. [51] to study the second and third cumulant in single and double quantum dot models. Explicit derivatives of (27) and the eigen-decomposition of the Kernel was used then to that end. In this subsection we show that only the stationary solution of the problem (solution to an algebraic equation) is needed to compute the finite-frequency current cumulants. We give analytical expressions (valid within the Markovian approximation) for the noise (second cumulant) and skewness (third cumulant) of the distribution of charge flowing through a conductor.
Let us decompose the Fourier transform in equation (18) into a set of Laplace transforms (defined as f (z) := ∞ 0 dte −zt f (t)), and the cumulant averages in terms of moments (c.f. for example Eq. (2.8) in Ref [62]). Doing this we find § The notation ">" denotes the unsymmetrized correlation function corresponding to the time ordering t N > . . . > t 2 > t 1 . Symmetrization in the frequency space implies adding the part corresponding to negative z and summing over all the possible switchings of frequencies. The subscript m means moment. These can be obtained as with z ≡ (z 1 , . . . , z N ) T and In the frequency domain, the prescription > can be taken similarly, that is z N > . . . > z 2 > z 1 , and finally symmetrize the result. and The advantage of having moment averages is that we can use a diagrammatic technique (see Appendix A) to easily obtain the correlation functions. Symmetrizing , we arrive to the following expressions for the current, noise and skewness of the current distribution (c.f. Appendix A for details): being Ω 0 (z) := [z − W(χ = 0)] −1 . These equations generalize the zero-frequency results found in Ref. [61] (c.f. their eqs Eqs. (7) and (8)) to finite frequencies (see Appendix A for the zero-frequency limit of (33)- (35)). Results for higher-order cumulants can be similarly obtained.
The relation between cumulants and moments can be formally expressed more generally at the level of the generating function. To do this one should follow the derivation by Kubo [62], making use of the property exp( i n i χ i ) = exp{ exp( i n i χ i ) − 1 c } in our context, arriving to a similar result to (7.25) in Ref. [62]. This allows for the calculation of frequency-dependent cumulants of the current distribution up to any order, reproducing in particular the results presented above. If a diagrammatic expansion in the Liouvillian space is used [68,69,70], cumulant averages become particularly useful since one can then keep only connected diagrams as those contributing to the average, in a similar way that this is done in quantum field theory.

FCS of total and accumulated currents
At finite frequencies, to have a theory consistent with current conservation it is essential to include the so called displacement currents [5]. This point is of vital importance to reproduce correctly the noise spectra measured experimentally. Although our discussion has focused, by construction, on particle currents so far, we show here how to include the effect of displacement currents in our formalism.
Let us illustrate this point by considering a quantum well with two terminals (L and R) in contact with Fermi leads at different chemical potentials. There will be then a net current flowing through both terminals, but also, charge can "accumulate" in the well for some time. Therefore charge conservation can be expressed aṡ with Q the charge in the well and I L , I R referring to the currents through the left and right contacts respectively.Q represents the displacement current, I dis , which can be partitioned as I dis = (α + β)I dis = I R dis + I L dis , where α and β describe how the displacement current is divided between the right and left reservoirs (obviously α + β = 1). This partitioning allows us to write the current conservation as In the simplest wide-band limit, the partition coefficients can be written in terms of tunnel rates only as α := Γ R /(Γ L +Γ R ) and β := Γ L /(Γ L + Γ R ) [71] + , and this will be the partitioning we will use throughout the paper.
Experimentally, one can measure correlations of the current through the device by transport measurements [20,21], or indirectly by studying the current through a charge sensor, such as a quantum point contact [22,23,24], that reveals whether the well is "charged" or "uncharged". The second method gives the statistics of the transport current only for very large bias voltages (unidirectional counting) but, in general the time-dependent transport current and the charge statistics are different. Morevover, when the device itself is used as a detector, the difference between transport fluctuations and charge fluctuations leads to profound physical consequences. Unlike the inelastic backaction induced by current fluctuations of the detector [65], the one induced by charge fluctuations is the fundamental Heisenberg backaction associated with the measurement [73]. Both transport and charge fluctuations can be accounted for in our formalism by considering counting fields which lead to respective jump operators where J χ,L and J χ,R refer to the two independent tunneling processes occurring at each barrier. The "total" cumulant through a two terminal device can be then calculated performing derivatives of the CGF with respect to χ tot defined in (37). This leads to expressions (33)-(35) with J ¶ In this paper we mean "total" cumulant when a subscript is ommited. + In a Coulomb Blockade model, the partition coefficients read α = C R /(C L + C R ) and β = C L /(C L + C R ), where C L , C R are the capacitances of each barrier and we have neglected capacitive effects from the gate. See for instance [72]. accum (ω). Notice that for a capacitive conductor, due to the relation between charge and voltage, this charge noise is proportional to the voltage noise. Finally the "left" and "right" cumulants can be computed with J 0 → J 0,L and J 0 → J 0,R respectively in (33)- (35).

Results
To illustrate our method, we analyze the transport statistics of the prototypical example of spinless electrons passing through a SRL model. The system consists on a two-terminal conductor with a discrete energy level in the central region, and is described by the Hamiltonian where k is the momentum and |0 and |1 are the only two possible states (referring to empty and occupied level) due to Coulomb blockade, with respective energies 0 and ε.
In the infinite bias limit (voltage V much larger than the other energy scales, excepting the bandwidth of the Fermi leads) the Hamiltonian (43) leads to the kernel expressed in the basis {|0 , |1 }, and where Γ α ≈ Γ α (E) := 2π k |V kα | 2 δ(E −ε kα ) are the rates accounting for the system-reservoir coupling. Using (33), (34) and (35), the simplicity of the model allows us to derive analytic results in this limit; for example the current gives I stat = Γ L Γ R /Γ, where Γ := Γ L +Γ R , and the "total" noise expressed in terms of the Fano factor (F (2) := S (2) /I stat ) reads At finite bias voltages, the kernel in Eq. (44) is no longer valid. Among the various choices to calculate W(χ) in this case, we use Schoeller's approach [68,69,70] (c.f. Appendix B) which allows us to calculate the Markovian kernel to the desired order (sequential tunneling in our case) without further uncontrolled approximations (such as the secular approximation). It is important to mention that the frequencydependent shot noise of the SRL model can be solved exactly [74], and therefore one does not need to use the above approximations. However, to the best of our knowledge, a finite-frequency study for this model beyond shot noise is yet lacking. Here we use the exact solution as a benchmark of the Markovian approximation in order to identify the regions of validity of our theory. This benchmark is important because most of the papers in the literature discussing shot noise in the context of QMEs make use of the Markovian approximation.
Another important check for the theory is to reproduce the fluctuation-dissipation theorem (FDT) in the appropriate regimes. Near linear response, that is, for applied voltages V much smaller than the temperature T , the low-frequency noise spectrum should follow the Jonhson-Nyquist relation S (2) = 2kT G [75,76], where G is the dc conductance. This equilibrium FDT was later extended by Callen and Welton to include quantum fluctuations [77], relevant when the measured frequencies are larger than the temperature. The FDT takes then the form S (2) (ω) = ωcoth( ω 2kT )G(ω), where G(ω) is the ac conductance. This expression can be equivalently written in terms of the Bose-Einstein distribution N BE (ω) ≡ 1/[e ω kT − 1], since coth( ω 2kT ) = 2N BE (ω) + 1 = N BE (ω) − N BE (−ω), and it becomes clear that the symmetrized noise, which we are considering here, contains both absorption and emission. Out of equilibrium, a fluctuation-dissipation relation can be also found for some particular cases, such as tunnel junctions [56,57] or for quantum dots in the weak cotunneling regime [58]. For a two-terminal conductor driven out of equilibrium, the symmetrized * noise spectrum takes the general form [78,79,32,5] where D n is the transmission coefficient of the conduction channel n. This expression has various interesting limits. First, in the tunneling regime (D n ≪ 1), it gives the non-equilibrium fluctuation-dissipation theorem (NEFDT) as reported in [56,57] for tunnel junctions: expression that is also exact for quantum dots in the weak cotunneling regime [58], and whose zero-frequency limit S (2) = I stat coth( eV 2kT ) has been derived in the context of counting statistics [55]. For low voltages, eV ≪ kT , equation (46) recovers the Callen and Welton equilibrium relation, and if also ω ≪ kT , it gives the Johnson-Nyquist FDT (thermal noise regime). Finally, if eV ≫ kT, ω (shot noise regime), we find S (2) = ζI, where the coefficient ζ ≡ n D n (1 − D n )/ n D n is the Fano-factor.
In this section we show the solution given by our theory for the SRL model, as well as how it recovers the the FDT and the NEFDT in the appropriate limits. By contrasting these results with the exact solution, we will be able to show that the Markovian approximation does not contain quantum fluctuations, thereby needing a non-Markovian approach to capture the physics of quantum noise [50]. It is therefore interesting to see to what extend the four results coincide, and in what regimes our Markovian approach is valid and recovers the proper physics. We will see that when kT ≫ eV, ω, the theory captures well both the exact solution and the FDT and NEFDT. Also, in transport configurations, with the level within the bias window, the Markovian approximation agrees well with the exact solution, reproducing in particular previous studies with eV ≫ kT [80,15,16,17]. However, in a situation with level outside the bias window, the Markovian approach presented here does not capture quantum noise physics, effect that we observe at high frequencies ( ω kT, eV ). Although in this situation transport due to cotunneling processes becomes more relevant, the difference is due to the Markovian assumption as can be seen using a non-Markovian extension of the theory [50]. We also study the finite-frequency skewness of the current distribution as given by equation (35). This shows to be insensitive to thermal fluctuations near equilibrium, therefore revealing * For a non-symmetrized version of the noise spectrum through a two-terminal conductor c.f. [65]. the "shot" contribution in a situation in which thermal fluctuations dominate in the noise spectrum (kT ≫ eV, ω).

Zero-frequency counting statistics.
Let us start by showing the zero-frequency noise spectra corresponding to the SRL model. Although this limit has already been studied in detail, a full comparison between our theory and the exact solution will help to understand the finite-frequency case. Particularly important is the linear response regime, at which studies of this model are scarce. As mentioned before, in this regime the noise should exhibit thermal fluctuations in order to fulfill a fluctuation-dissipation relation, while the skewness, on the other hand, should go to zero as the voltage V goes to zero [55]. Fig. 1 shows how our calculation captures correctly the FDT, S (2) = 2kT G, in the proximity of linear response, kT eV . For comparison, we also plot the zero-frequency limit of the NEFDT in Eq. (47), namely S (2) = I stat coth( eV 2kT ) [55]. In the opposite regime, kT eV , the Markovian approximation is larger than the exact solution, discrepancy that can be understood as originated from the lack of cotunneling contributions in our calculation [70]. As expected, below kT /eV ∼ 1 the FDT is not fulfilled. We can also see that the NEFDT, exact for tunnel junctions, for a two-terminal device performs quite badly when kT eV , but correctly in the opposite limit. This failure of the NEFDT at low temperatures disappears when the level is outside the bias window. This is a low-current regime, and thus a tunneling limit. The inset of Fig. 1  range of temperatures. At finite frequencies we will expect a quantum noise step in the spectrum at frequencies ω ∼ ε, effect that will be studied in the next subsection.
The behaviour of the zero-frequency noise spectrum close to equilibrium with respect to ε is shown in Fig. 2a. Here we see how the FDT is fulfilled by our theory and the good agreement with the exact solution. We also plot the zero-frequency skewness, that although of small magnitude in the same scale, is nonzero in a situation where the noise spectrum is completely dominated by thermal fluctuations. This insensitivity of the skewness to temperature allows us to extract intrinsic correlation effects in near-equilibrium conditions. In Fig. 2b we show the same quantities as a function of voltage. Interestingly, the Markovian result coincides with the exact solution in the whole range of voltages. The FDT, however, starts to disagree with these for voltages eV /ε 0.2. As anticipated, the skewness vanishes as the voltage goes to zero. In Fig. 3a we plot noise and skewness as a function of ε for increasing voltages. As the bias increases, the skewness (dashed lines) shows peaks evolving into plateaus at values of ε corresponding to the chemical potentials of the reservoirs. This effect, which is due to non-equilibrium fluctuations, is completely masked in the noise (solid  lines) even at the highest voltages due to thermal fluctuations. This is clearly seen in Fig. 3b, where we show the same comparison after substracting thermal fluctuations to the noise value (excess noise defined as S (2) (V ) − S (2) (V = 0)). Here it is clear that at low detuning, ε eV /2 (position of the peaks in the figure), and when kT eV , the skewness can reveal the "shot" contribution, while this is masked by thermal fluctuations in the noise spectrum.

Finite-frequency counting statistics.
To study the case of finite frequencies, we use our formulae (33)- (35) applied to the SRL model. In a situation with the level within the bias window, we find a similar behaviour to Fig. 1. However, now the NEFDT -equation (47) -is fulfilled for temperatures kT eV + ω. This is shown in Fig. 4. Remarkably, at finite frequencies the Markovian approximation is basically exact in this direct transport regime. In the high-bias regime eV ≫ ω, kT we also find that the Markovian approximation agrees perfectly with the exact result (not shown), in accordance with previous studies [80,15,16,17]. When the level lies outside the voltage window, the situation changes drastically (see inset of Fig. 4). Here the Markovian approximation is no longer appropriate when kT eV + ω. Both the exact result and the NEFDT contain quantum fluctuations, while the Markovian calculation only captures the thermal contribution (and therefore fulfills the FDT). The exact solution presents a small structure at temperatures of the order of ε. This cannot be resolved with the NEFDT. As expected, when kT eV + ω, all curves coincide.
The general trends explained so far become more evident in Fig. 5, where we plot  In the inset all fluctuations are thermal and, therefore, the shot noise and the FDT coincide in the whole range of temperatures. When kT ≤ eV + ω, the NEFDT is above due to quantum fluctuations. The exact solution contains also corrections due to cotunneling processes, which are dominating in this regime.
noise and the fluctuation-dissipation theorem near linear response as a function of frequency. Here the Markovian noise is always flat and equals S (2) = 2kT G, whereas the NEFDT and the exact solution lie above and show quantum noise steps. Let us start by considering the case ε = 0. In the whole range of frequencies, the Markovian approximation is basically exact in this situation of direct transport. The NEFDT shows the correct zero-frequency limit, since then the fluctuations are purely thermal. At high frequencies, however, the NEFDT converges to the Poisson value of a single barrier with tunneling rateΓ/2, beingΓ := (Γ L + Γ R )/2, (c.f. S (2) L in the plot). This is in agreement with the validity of equation (47) for a tunnel junction, and in contrast with the exact solution, which contains partitioning, and therefore its ω → ∞ limit isΓ/4. In the case in which the energy lies outside the bias window (ε/kT = 5 in the figure), transport is possible because of the finite temperature as well as due to quantum fluctuations. The Markovian noise only contains the former and is flat with frequency, while the exact result contains both and shows a quantum noise step centered at ω/kT = ε/kT = 5. Although in this regime cotunneling contributions are important, the difference lies in the Markovian approximation, as can be seen using a non-Markovian theory [50]. The NEFDT shows in this situation a quantum noise step centered at ω = 2ε. This discrepancy with respect to the exact solution can be understood in terms of the tunneling approximation leading to Eq. (47), which presents a step located at an effective chemical potential 2ε. Again, the high frequency limit coincides with that of S  47)) and the exact solution. S (2) (ω) is flat for the whole range of frequencies, and coincides with the equilibrium FDT as expected. The NEFDT however disagrees with these two, showing also quantum fluctuations which are absent in the Markovian noise spectrum. The quantum noise steps shown by the NEFDT are however at ω = 2ε, in contrast to the exact solution, which shows steps at ω = ε. This is due to the fact that the NEFDT works well for tunnel junctions, but does not capture partition noise. This becomes clear also from the saturation value at large frequencies, as described in the text. Rest of parameters: hΓ L /kT = hΓ R /kT = 0.05. The inset compares the exact solution with the Markovian approximation for a different regime, namely eV /kT = 25. We see that while the Markovian limit is flat for all frequencies, the exact solution presents a dip at ω = ±|ε ± eV /2|. Rest of parameters: ε = 0, hΓ L /kT = hΓ R /kT = 0. 25. solution, whereas the exact solution presents a dip at ω = ±|ε ± eV /2| (coinciding with the position of the chemical potentials with respect to the energy level). This clearly illustrates how the Markovian approximation captures well the physics in the linear response regime and a direct transport configuration, but when the frequency is comparable to the applied bias, it fails to capture the quantum noise.
We now proceed to discuss the finite-frequency noise and skewness spectra of the total-and particle-current distribution. In the previous discussion, the Markovian noise was always flat as a function of frequency, something which is well known for symmetric systems (Γ L = Γ R ) -see for instance [5]. In the case Γ L = Γ R , a proper partitioning of displacement currents (see discussion in subsection 2.5) becomes essential as we will show next, and the way this is made affects significantly the spectrum. Fig. 6a shows S (2) (ω) in a transport configuration, eV /hΓ R = 50 and ε = 0. As in the results shown previously, the total noise spectrum is flat for a symmetric configuration. Interestingly, this flat behaviour persists even when the system is made asymmetric (Γ L = Γ R ). This is due to the current-partitioning model assumed here: I tot = αI L + βI R with α = Γ R /Γ and β = Γ L /Γ. The noise spectrum corresponding to particle currents displays information about the rates; in contrast to the total noise, it shows a dip with half-width 2Γ. In Fig. 6b we show the skewness along the representative direction ω ′ = −ω. Interestingly, the skewness corresponding to the total current starts to develop a dip that shows the asymmetry of the system. The particle-current skewness presents a similar behaviour to the noise counterpart.  Figure 6. a) S (2) (ω) and b) S (3) (ω, −ω) for eV /hΓ R = 50, kT /hΓ R = 20, and ε = 0. The spectra for particle currents and for total currents significantly deviate form each other, even for large asymmetry. The insets correspond to noise and skewness through the left barrier for Γ L /Γ R = 10.
However, for (3− √ 5)/2 Γ L /Γ R (3+ √ 5)/2 it develops a minimum whose position depends on the value of the rates. In the asymmetric case, Γ L = Γ R , the particlecurrent noise and skewness can even present different lineshapes. This can be seen contrasting the insets of Fig. 6a and Fig. 6b. In the linear response regime the curves for the noise look similar (not shown). The skewness on the contrary goes to zero in magnitude and shows a structure that depends on temperature, and that changes from a dip to a peak as ε is increased from zero to a finite value. In summary, we see that the spectra for total and particle currents differ significantly from each other even for large asymmetry. This means that the assumption of calculating noise spectra using particle currents only, used commonly in the literature, is flawed. Here we have assumed the current partitioning given by α = Γ R /Γ, β = Γ L /Γ. If the more simplistic partitioning α = β = 1/2 is assumed, the results for the total cumulants in the asymmetric case change drastically (not shown). In particular, the noise is then no longer flat but has a dip structure, and the skewness shows a peak around zero frequency.

Conclusions
In conclusion, we have developed a theory of frequency-dependent counting statistics of electron transport through nanostructures within the framework of Markovian quantum master equations. We have illustrated our method with calculations of noise and skewness in a single resonant level model at finite bias voltages and frequencies.
By comparing with both the exact solution and the finite-frequency version of the nonequilibrium fluctuation-dissipation theorem, Eq. (47), we have identified the regimes of validity of our Markovian theory at finite frequencies. In particular, we have shown that the Markovian limit is basically exact in transport configurations (level within the bias voltage window), as long as ω ≫ eV or ω ≪ eV . In intermediate situations where ω ∼ eV , or with the level outside the bias window, the Markovian limit fails at finite frequencies due to the lack of quantum fluctuations [50].
We have also discussed how the noise spectra for particle currents and for total currents significantly deviate from each other, even for large asymmetries Γ R /Γ L = 1. This demonstrates that calculating spectra using particle currents only leads to incorrect results in general. Our method allows the calculation of finite-frequency current cumulants of arbitrary order, as we have explicitly shown for the second and third order cumulants, Eqs. (34) and (35). These formulae generalize previous zerofrequency expressions and can be viewed as an extension of MacDonald's formula beyond shot noise. Recently, this has been extended to study the time-averaged shot noise spectrum in the presence of periodic ac fields [81,82]. Interesting extensions of our study along these lines would allow us to study frequency-dependent highorder cumulants in nanostructures driven by time-dependent fields, or, even more challenging, in systems showing nontrivial non-linear dynamics such as self-sustained oscillations without external time-dependent driving [83].
0 Ω 0 (z 23 )Ω 0 (z 123 )J (1) 0 where z ij := z i + z j , z ijk := z i + z j + z k . Next we use and add the "lesser" part (<) corresponding to negative Laplace frequencies. At this point we change to "physical" frequencies ω := ω 2 + . . . + ω N , ω ′ := ω 3 + . . . + ω N , etc., and symmetrize the result. This means adding the expressions corresponding to all the possible frequency switchings. In this step we take into account that where η → 0 is a small parameter coming from the "greater" (>) or "lesser" (<) parts. Finally, we make use of this energy conservation inherited from the time-translational symmetry of the cumulants. We then arrive to the equations (33)-(35) used in the main text. Importantly, after frequency symmetrization one can realize that the first three cumulant formulae are equal to their moment counterparts.
Interestingly, the results (A.1)-(A.3) given above can be derived following a diagrammatic technique, similarly to how this is done with Feynman diagrams in the expansion of the partition function or the S-matrix. Similarly we write the CGF in terms of a series expansion, either in the time domain or in the frequency space.
To that end we expand each of the χ-dependent propagators in the CGF as a Dyson series: This suggests the use of diagrams of the form given in Fig. A1 (a). In the frequency domain ♯ this rules are: ♯ If we work in the time domain, it is enough to label the propagating lines with the corresponding times at the beginning and end of each line (see Fig. A1 (a)).
• To each bare propagator Ω 0 (z k ) in the expansion we associate a line with a superscriptk ≡ N i=N +1−k i, where N is the order of the cumulant we want to obtain.
• To each jump operator Jχ k in the expansion we associate an encircled cross with superscriptk. The formula for the generating function to a given order can therefore be written diagrammatically. For example, to second order we have (A.7) We can then multiply the different terms using diagrams as described above. The multiplication of propagators implies joining them together. The result can be simplified using the property J χ1+χ2 = J χ1 J χ2 + J χ1 + J χ2 , which diagrammatically is denoted as Here, the super-index 12 denotes an associated frequency z 1 + z 2 and counting field χ 1 + χ 2 . Next, to arrive to the frequency-dependent moment, we take derivatives with respect to counting fields. Diagrammatically, the derivative ∂ χ k means removing a circle with index k. From here we can rewrite the expression analytically. The outcome can be simplified using (A.4), and needs to be multiplied by z 1 z 2 (case N = 2), coming from the Fourier transform of the time derivatives in the frequency domain. We finally need to take the average in the stationary state and symmetrize the result as dictated by T S . With the diagrammatic approach we realize that the diagrams contributing to the final result can be arranged in tables (see Fig. A1 (b) to (d)). These reproduce the results given in (A.1)-(A.3). To construct these tables we proceed according to the following rules: • To arrive to an expression for the cumulant of order N , write a table with N time (frequency) intervals and corresponding superscriptsk. The propagation of time will be taken from right to left. • Write all the possible diagrams having N crosses (jumps) distributed in the different intervals, with the constraint that the maximum number of crosses in each is set by the corresponding indexk. Diagrams with n jumps occurring at the same time have to be included as well. These crosses are enclosed together with a box, and contribute with the jump operator J (n) 0 := ∂ n χ J χ | χ=0 . • Taking into account that jumps occurring in the same interval are indistinguishable, and that each cross can be associated to one of the possible counting fields χ 1 , . . . , χ k present in that interval, write the multiplicity of each diagram on the right. • Write the mathematical expression corresponding to each diagram (see Fig. A1 (a)) and sum the different terms evaluated at z = iω. • Take the average in the stationary state, and multiply by (−i) N . The resulting expression corresponds to the unsymmetrized ("greater", >) moment of the number of particles.
• Finally, symmetrize the result, adding all the possible switchings of frequencies.
This gives the symmetrized moment of the current distribution. The result can be simplified using (A.4) and (A.5).
As mentioned above, explicit derivation gives the same result for the expressions of cumulants of the current distribution as those derived for the moments up to N = 3. To higher orders it is unknown for us if this property still holds or not. Expressions for the cross correlations, e.g. S (2) LR := I(t 1 )I(t 2 ) c , between two (or more) stochastic processes, e.g. L and R, can also be derived with this technique. To this end we simply need to label each of the jumps occurring at L or R accordingly (see Fig. A1 (c)), having two types of jump operators, J L and J R . Also, expressions for the total current (αI L + βI R ) or accumulated current (I L − I R ) can be derived using the jump operators (39)- (40) in the diagrams.

A.2. Equivalent form.
We can write down an equivalent form to expressions (33)- (35). This will allow us to obtain an analytical expression for their zero-frequency limit, which is not well defined in the form given above. To this end we make use of the projectors P := |0 0| and Q := 1 − P , where P projects onto the subspace spanned by the stationary state |0 ≡ ρ stat S † †; and we define the pseudo-inverse R 0 (z) := QΩ 0 (z)Q, such that Ω 0 (z) = R 0 (z) + P/z.
(1) 0 † † The state 0| denotes the left eigenvector of the Liouvillian. The tilde indicates that it is not the adjoint to |0 , since the Liouvillian is not Hermitian. (A.11) Now we make use of the delta function to write z 2 = −z 1 in the noise expression and z 3 = −z 1 − z 2 in the skewness result. Performing the change of variables z 1 → −iω, z 2 → iω in the noise and z 1 → −iω, z 2 → iω − iω ′ , z 3 → iω ′ in the skewness, we obtain iI stat = J (A.14) The limit ω → 0 of these expressions is well defined, and they can therefore be used to check that the proper result is recovered in that limit.
As mentioned, expressions (A.12)-(A.14) are well behaved when ω → 0. The zerofrequency noise comes straightforwardly from (A.13) setting ω = 0. For the skewness, this limit requires nevertheless noticing that Which is the zero-frequency limit found in [61].  Figure A1. Diagrammatics to obtain the frequency-dependent cumulants. a) Building pieces for the diagrammatic technique. A line is associated to a bare propagator, a cross with a jump operator, and a circle with the time dependence of J (term (e iχ − 1) in the single-particle unidirectional tunneling case). Derivatives with respect to counting fields eliminate circles correspondingly. Diagrams can be simplified using rules like (A.8). b) Diagrams for the noise.