Non-renewal statistics in quantum transport through the eyes of first-passage and waiting time distributions
Samuel L. Rudge Daniel S. Kosov




College of Science and Engineering, James Cook University, Townsville, QLD, 4811, Australia


The waiting time distribution has, in recent years, proven to be a useful statistical tool for characterising transport in nanoscale quantum transport. In particular, as opposed to moments of the distribution of transferred charge, which have historically been calculated in the long-time limit, waiting times are able to detect non-renewal behaviour in mesoscopic systems. They have failed, however, to correctly incorporate backtunneling events. Recently, a method has been developed that can describe unidirectional and bidirectional transport on an equal footing: the distribution of first-passage times. Rather than the time between successive electron tunnelings, the first-passage refers to the first time the number of extra electrons in the drain reaches +1. Here, we demonstrate the differences between first-passage time statistics and waiting time statistics in transport scenarios where the waiting time either cannot correctly reproduce the higher order current cumulants or cannot be calculated at all. To this end, we examine electron transport through a molecule coupled to two macroscopic metal electrodes. We model the molecule with strong electron-electron and electron-phonon interactions in three regimes: (i) sequential tunneling and cotunneling for a finite bias voltage through the Anderson model, (ii) sequential tunneling with no temperature gradient and a bias voltage through the Holstein model, and (iii) sequential tunneling at zero bias voltage and a temperature gradient through the Holstein model. We show that, for each transport scenario, backtunneling events play a significant role; consequently, the waiting time statistics do not correctly predict the renewal and non-renewal behaviour, whereas the first-passage time distribution does.


Non-renewal statistics in quantum transport through the eyes of first-passage and waiting time distributions
    Daniel S. Kosov
    December 30, 2023
============================================================================================================



§ INTRODUCTION

Molecular electronic devices generally consist of some quantum system, such as a molecule or quantum dot, sandwiched between two macroscopic metal electrodes <cit.>. One of the main objectives of mesoscopic physics is to use the particle transport through such devices to understand the nature of these quantum systems. Directly measuring transitions internal to the quantum system destroys any quantum coherence it may exhibit; instead, one must reconstruct the system behaviour from transport properties available from measurements of bath variables, such as the stationary particle current ⟨ I ⟩. Molecular devices are small enough, however, that thermal effects and stochastic single electron transport cause fluctuations in ⟨ I ⟩. Such fluctuations hold a dual interest in mesoscopic electron transport. On the one hand they impede the performance of electronic components that require reliable currents; but on the other hand they are largely determined by the molecule's underlying dynamics and can thus be used to probe quantum behaviour, which is difficult to observe. Consequently, multiple theoretical tools describing fluctutations in quantum electron transport have been developed <cit.>. 

The full counting statistics (FCS) calculates cumulants ⟨⟨ I^k⟩⟩ of the distribution of total transferred charge n in a fixed time interval. <cit.> The second current cumulant, the zero-frequency noise, can reveal the internal dynamics of the molecule <cit.>. For example, super-Poissonian shot noise has many origins: it can arise from telegraphic switching due to spin-polarised leads <cit.> or inelastic cotunneling <cit.>; negative differential resistance due to asymmetric coupling <cit.>; the dynamical channel blockade <cit.>; or avalanching electrons due to interactions with a vibrational mode <cit.>. Higher order cumulants are necessary to fully characterise the transport when the current distribution is non-Gaussian <cit.>. The major limitation of fixed time tools like the FCS is that they have historically been calculated in the infinite time, or zero-frequency, limit <cit.>. Although finite frequency, and thus finite time, FCS have been developed <cit.>, it is numerically and analytically easier to work in the zero-frequency regime. In recent years, however, statistical tools that operate in fluctuating time but fixed n have arisen, motivated by the long time limitations of FCS.

The most prominent flucutating time tool is the waiting time distribution (WTD) w(τ): the conditional probability density that, given an electron has tunneled to the drain electrode, the next electron will tunnel to the drain after a time interval τ <cit.>. Since the inception of waiting times in mesoscopic electron transport, a common question has been whether the WTD contain information distinct from the FCS. It has been shown that the large deviation rate functions describing zero-frequency FCS are equivalent to those for waiting times <cit.>; the WTD cumulants ⟨⟨τ⟩⟩ can be combined to exactly reproduce the zero-frequency current cumulants. <cit.> These relations only hold, however, when the renewal assumption is satisfied: w(τ_1,τ_2)=w(τ_1)w(τ_2), where w(τ_1,τ_2) is the joint probability for successive waiting times. When subsequent waiting times are correlated the renewal assumption is violated and these relationships no longer hold; the WTD contains information separate to the FCS. Interestingly, even in Markovian systems non-renewal behaviour can arise from telegraphic switching <cit.>, superconductivity <cit.>, quantum coherence <cit.>, and the presence of strong interactions <cit.>. WTD measurements are experimentally accessible via sophisticated real-time single-electron detection techniques <cit.>; or alternatively the WTD may be reconstructed from low-order current correlation functions <cit.>.

The WTD in quantum electron transport can be calculated via scattering theory <cit.> and non-equilibrium Green's functions <cit.>. In order to connect it with other fluctuating time statistics, however, we focus on WTDs calculated from quantum master equations, which rely on quantum jump operators defined from the Liouvillian <cit.>. Unfortunately, WTDs calculated in this manner are so far only able to reproduce ⟨⟨ I^k⟩⟩ (k>1), as according to renewal theory, for unidirectional transport. For bidirectional transport, one may define WTDs between tunnelings to and from the drain; the total average current, regardless of whether the transport is renewal, is then ⟨ I ⟩ = 1/⟨τ⟩_F - 1/⟨τ⟩_B, where ⟨τ⟩_F is the average waiting time between forward tunnelings, conditioned upon the initial probability of a tunneling to the drain, and ⟨τ⟩_B is the corresponding backward tunneling average waiting time. It is not clear, however, how to obtain the higher order total current cumulants from these two WTDs. In some cases, furthermore, the processes that contribute to the current do not appear in the Liouvillian; for example, elastic cotunneling transfers electrons from the source to the drain in the same quantum process without altering the state of the intermediate system <cit.>. To fully describe the transport in these cases, one must use the n-resolved master equation <cit.>; but waiting times calculated from this method <cit.> are only defined for unidirectional transport <cit.>. Bidirectional transitions play an important role in mesoscopic electron transport <cit.>, so it is desirable to have a fluctuating time statistic applicable to this regime.

Consider instead the first-passage time distribution (FPTD): the distribution of times until the net number of transferred electrons in the drain first reaches some number n. First-passage times have recently gained traction in non-equilibrium thermodynamics; first-passage time theory has been developed for fluctuations in entropic variables describing stationary Markovian processes <cit.>. This work, and subsequent experimental and theoretical research from Singh et al. <cit.>, has predicted and verified fluctutation relations for first-passage times. Saito and Dhar <cit.> have developed first-passage time notation for Markovian quantum systems, which is suitable for mesoscopic electron transport described by n-resolved master equations. Ptaszynski <cit.> has also recently brought first-passage times into line with FCS and WTDs, by developing a method to calculate correlations between subsequent first-passage times; deriving exact relationships between the FPTD cumulants and the FCS cumulants; and demonstrating that these relationships disappear when the renewal assumption is violated. The expectation now is that first-passage times are able to analyse the non-renewal behaviour of molecular systems similarly to how WTDs have been used, except that first-passage times are applicable to bidirectional transport. Considering these developments, it is interesting to compare the WTD and the FPTD in systems and regimes where one expects bidirectional transitions. 

This paper's purpose is to demonstrate the usefulness of FPTDs by comparing them to the WTD in regimes where the WTD either cannot be calculated or does not correctly recreate the FCS. We demonstrate this in three transport scenarios. First, cotunneling through an Anderson impurity is considered. Several authors have theoretically <cit.> and experimentally <cit.> shown that for a particular voltage range super-Poissonian noise arises from telegraphic switching induced by inelastic cotunneling processes. In this voltage range backtunneling processes can occur; so until now it has remained impenetrable to waiting times. One would expect that telegraphic switching would be accompanied by positive correlations in the subsequent waiting times, as the transport switches between channels with different characteristic ⟨τ⟩; however, we demonstrate that the transport in this regime is actually characterised by small negative correlations in the first-passage times. Second, we examine sequential tunneling through the Holstein model due to a voltage bias. Waiting times have already been theoretically calculated for this scenario <cit.>; we demonstrate that those calculations do not correctly reproduce the FCS and that the correlations between the first-passage times are slightly smaller than those between the waiting times. Finally, we consider sequential tunneling through the Holstein model with zero bias voltage but with a temperature gradient: a scenario in which backtunneling events majorly contribute to the transport <cit.>.

The paper is organised as follows. Section <ref> outlines the three relevant statistical tools: FCS, the WTD, and the FPTD. Section <ref> describes the dynamics of the three transport scenarios, written in the form of Markovian master equations. Section <ref> contains results from the statistical tools calculated for the three transport scenarios, while Section <ref> summarises the work and the Appendix details cotunneling calculations and vibration jump operators.

Throughout the manuscript we use natural units: ħ=k_B=e=1.





§ STATISTICAL TOOLS
 



 §.§ n-resolved master equation
 

We consider a general quantum system weakly coupled to a number of thermal reservoirs, which is described by a Markovian master equation for the reduced system density matrix:


    𝐏̇(t)     = 𝐋𝐏(t),

where {𝐏}_l(t) is the probability that the system is in discrete state l and 𝐋 is the Liouvillian containing the system dynamics. Throughout the paper, we refer to this as the standard master equation, to distinguish it from the n-resolved master equation. We ignore coherence between system states, so that Eq.(<ref>) reduces to a rate equation written in the system basis. The system moves from state l to state k via tunneling of particles or quasi-particles from or to the thermal reservoirs, and the transition rate for this process is Γ_kl. The Liouvillian off-diagonals are [𝐋]_kl = Γ_kl and the diagonals are [𝐋]_kk = -∑_l≠ kΓ_lk. 

The system starts in some state 𝐏(0), so that the general solution of Eq.(<ref>) is


    𝐏(t) = e^𝐋t𝐏(0).

We assume that the system has a unique stationary state 𝐏̅, which is the null vector satisfying 𝐋𝐏̅ = 0.

In order to calculate the FCS, FPTD, and WTD, it is necessary to resolve the master equation on the net number of particles n transferred to the drain electrode:


    𝐏̇(n,t)     = ∑_n'𝐋(n-n')𝐏(n,t),

where {𝐏(n,t)}_k now represents the probability that, at time t, the net number of extra electrons in the drain is n and the system is in state k. It is important to note that despite the term “extra electrons” n may be zero, a negative integer, or a positive integer. From here, following Ptaszynski's notation <cit.>, we will refer to n as the jump number. Sequential and cotunneling processes add or remove at most one electron to or from the drain, so n-n' = 0,± 1. The Liouvillian can then be split into two jump operators containing all the transitions that move electrons forward to the drain 𝐉_F and all the transitions that move electrons back from the drain 𝐉_B, and a part containing the remaining dynamics 𝐋_0 = 𝐋 - 𝐉_F - 𝐉_B. The n-resolved master equation is then


    𝐏̇(n,t)     = 𝐋_0𝐏(n,t) + 𝐉_F𝐏(n-1,t) + 𝐉_B𝐏(n+1,t).
 
The Fourier transform of 𝐏(n,t) 


    𝐏(χ,t)     = ∑_ne^inχ𝐏(n,t), and conversely
    𝐏(n,t)     = 1/2π∫_0^2π dχ e^-inχ𝐏(χ,t),

conveniently introduces a counting field χ that transforms Eq.(<ref>) from an infinite set of coupled differential equations to a solvable problem <cit.>:


    𝐏̇(χ,t)     = 𝐋(χ)𝐏(χ,t),  where 
    𝐋(χ)     = 𝐋_0 + 𝐉_Fe^iχ + 𝐉_Be^-iχ.


Now, 𝐏(χ,t) time-evolves according to the solution of Eq.(<ref>), with the initial condition 𝐏(χ,0) = 𝐏̅ since it is assumed that counting begins at t=0 when the system has already reached the steady state:


    𝐏(χ,t)     = e^𝐋(χ)t𝐏̅.




 §.§ Full counting statistics
 

The generating function for the cumulants of transferred charge is 


    M(χ,t)     = ln∑_n P(n,t) e^inχ,
 
where P(n,t) = [𝐏(n,t)] is the probability that n extra electrons are transferred to the drain in the time interval [0,t]. Throughout the paper we use the superoperator notation [𝐀] = (𝐈,𝐀), where 𝐀 is a column vector of length m and 𝐈 = (1,1,1,...) is a row vector of length m. Consequently, Eq.(<ref>) transforms to 


    M(χ,t)     = ln∑_n(𝐈,𝐏(n,t)) e^inχ,
 
with the inverse Fourier transform 


    (𝐈,𝐏(n,t))     = 1/2π∫_0^2π dχ e^-inχ e^M(χ,t).
 

Examining the trace of Eq.(<ref>), one sees that the part corresponding to the cumulant generating function is 


    M(χ,t)     = ln(𝐈,e^𝐋(χ)t𝐏̅).
 

Theoretically, one can now calculate all current cumulants in any time interval by direct differentiation of M(χ,t) in the limit χ→0: ⟨⟨ I^k⟩⟩ = 1/t(-i)^k∂^k/∂χ^kM(χ,t)|_χ=0, where t is the measurement time. The second cumulant, for example, is ⟨ I^2⟩ - ⟨ I⟩^2. It is well known that, except for the simplest systems, these calculations remain too difficult at finite t. Bagrets and Nazarov have shown that in the long-time limit, the cumulant generating function reduces to 


    M(χ,t)     = tΛ_max,

where Λ_max is the eigenvalue of 𝐋(χ) with the largest real part  <cit.>.

We focus on the Fano factor, the second cumulant scaled by the average current, which measures the zero-frequency noise power 𝒮(0) compared to that of a Poissonian process:


    F     = ⟨⟨ I^2⟩⟩/⟨ I ⟩
        = 𝒮(0)/2⟨ I ⟩.


For Poissonian processes F=1; sub-Poissonian and super-Poissonian processes are characterised by F<1 and F>1, respectively. 



 §.§ Waiting time distribution
 

If instead one wishes to examine the distribution of waiting times between subsequent tunnelings to the drain, then one may use the formula developed by Brandes <cit.>:


    w_F(τ)     = (𝐈,𝐉_Fe^(𝐋-𝐉_F) τ𝐉_F𝐏̅)/(𝐈,𝐉_F𝐏̅),

The waiting time distribution between tunnelings from the drain, w_B(τ), is similarly defined. 

There are tunneling events that do not appear in the standard master equation, yet change the drain number n. To include these events in the WTD, one must use the n-resolved master equation and the definition of the WTD from the idle time probability <cit.>:


    w_F(τ)     = 1/p∂^2/∂τ^2Π(τ),

where p is the probability for observing a tunneling to the drain and Π(τ) is the idle time probability: the probability that no tunneling events to the drain are observed by time τ. 

In the case of unidirectional transport, when 𝐉_B = 0 and all backtunneling transitions are safely neglected from the Liouvillian, p = ⟨ I ⟩ and the idle time probability is obtained from n-resolved master equation:


    Π(τ)     = lim_iχ→ -∞(𝐈,e^𝐋(χ)τ𝐏̅).

Now, 𝐋(χ) = 𝐋_0 + 𝐉_Fe^iχ and no backtunneling transition rates are included in 𝐋_0. Evidently, backtunneling events will be accompanied by an e^-iχ factor that will diverge in the limit iχ→ -∞. Of course, if the transport is unidirectional in the backward direction as opposed to the forward direction, then the same definition may be used with the limit iχ→∞; either way only one direction of tunneling events may be included. Using Eq.(<ref>) in the waiting time definition from Eq.(<ref>) yields the same form of the WTD as developed by Brandes, except now the exponent 𝐋(0) - 𝐉_F is obtained from the n-resolved master equation; when the transport is unidirectional the two methods coincide. The only difference between them is that when the jump operators can be defined in terms of a standard master equation Brandes' method can calculate both w_F(τ) and w_B(τ), whereas Eq.(<ref>) cannot.

The Laplace transform of the WTD is


    w̃(z)     =  (𝐈,𝐉(z - 𝐋_0)^-1𝐉𝐏̅)/(𝐈,𝐉𝐏̅),

from which the cumulants are calculated:

    ⟨⟨τ^k⟩⟩    = (-1)^kd^k/dz^klnw̃(z)|_z=0.

We have included the notation 𝐋_0 = 𝐋 - 𝐉 as shorthand for the exponent in Eq.(<ref>). 

The average current can be reconstructed from the average waiting times:


    ⟨ I ⟩ = 1/⟨τ⟩_F - 1/⟨τ⟩_B,

where ⟨τ⟩_F and ⟨τ⟩_B are the average forward and backward waiting times, respectively. Waiting times calculated from an n-resolved only satisfy Eq.(<ref>) for unidirectional transport<cit.>. The equality in Eq.(<ref>) holds regardless of whether the transport is renewal, which follows from examining the forward and backward currents separately. The forward current, reconstructed from the WTD, is


    ⟨ I ⟩_F    = ∑_k=1^∞k/⟨τ_1+τ_2++τ_k⟩_FP(k)

where P(k) is the probability that the forward current is k/⟨τ_1+τ_2++τ_k⟩_F. The average waiting times are defined via the joint probability distribution:



    ⟨τ_1+τ_2++τ_k⟩    = ∫ dτ_k∫ dτ_2∫ dτ_1(τ_1 + τ_2 + τ_k) w(τ_1,τ_2,,τ_k) 
        = ∑_p=1^k∫ dτ_p τ_p w(τ_p) 
        = k⟨τ⟩,


where we have made no assumption on the renewal nature of the transport. Eq.(<ref>) then reduces to 

    ⟨ I ⟩_F    = 1/⟨τ⟩_F,

assuming that {P(k)} spans the entire current space. The same reasoning follows for the backward current. 

The higher order current cumulants, in comparison, are only reconstructable from the WTD for unidirectional transport; for example,  the randomness parameter 


    R     = ⟨⟨τ^2⟩⟩/⟨τ⟩^2,


is equal to the Fano factor defined in Eq.(<ref>) for unidirectional transort satisfying the renewal assumption. The Pearson correlation coefficient measures the strength of correlations between successive waiting times:


    p     = ⟨τ_1τ_2⟩ - ⟨τ⟩^2/⟨⟨τ^2⟩⟩.


When p is nonzero, the renewal assumption is violated and the equality between Eq.(<ref>) and Eq.(<ref>) no longer applies. 

The first moment of the joint distribution ⟨τ_1τ_2⟩ is calculated similarly to the cumulants:

    ⟨τ_1τ_2⟩    = ∂/∂ z_1∂/∂ z_2w̃(z_1,z_2)|_z_1=z_2=0
        =  (𝐈,𝐉𝐋_0^-2𝐉𝐋_0^-2𝐉𝐏̅)/(𝐈,𝐉𝐏̅).




 §.§ First-passage time distribution
 

Unlike the WTD, which calculates the time interval between successive electron tunnelings to the drain, the FPTD calculates the time interval between a tunneling to the drain and when the jump number next reaches a specified integer n. Necessarily, the FPTD includes backward tunneling events, and it is defined from the n-resolved master equation. To derive the FPTD we follow the method neatly encapsulated by Saito and Dhar in the framework of Markovian master equations <cit.>, and explicitly explained by Ptaszynski <cit.>.

We return to the χ-dependent Liouvillian defined in Eq.(<ref>) and start the time-evolution in some arbitrary state 𝐏(0). The n-resolved probability vector is written using Eq.(<ref>) and Eq.(<ref>):


    𝐏(n,t)     = 1/2π∫_0^2π dχ e^-inχ e^𝐋(χ)t𝐏(0) 
        = 𝐓(n|t)𝐏(0).


Now, 𝐓(n|t) = 1/2π∫_0^2π dχ e^-inχ e^𝐋(χ)t is the evolution operator for 𝐏(n,t), which in Laplace space is 


    𝐓̃(n|z)     = 1/2π∫_0^2π dχ e^-inχ[z - 𝐋]^-1.

The matrix element [𝐓(n,t)]_kl is the conditional probability that, given the system is initially in state l at time t=0, it will be in state k at time t and the jump number will have reached n for the first time. Saito and Dhar <cit.>, Ptaszynski<cit.>, and Singh et al.<cit.> have all found exact expressions for 𝐓̃(n|z) for simple single reset systems, in which a tunneling to the drain leaves the molecule empty. In general, though, evaluating Eq.(<ref>) requires a numerical contour integration, which is the approach we take. 

Alternatively, 𝐏(n,t) is obtainable by considering the column vector of first-passage time probabilities 𝐅(n|τ). Here, [𝐅(n|τ)]_k is the probability density that at time τ the system is in state k and n extra electrons have been recorded for the first time in the drain. Premultiplication by the evolution operator 𝐓(0,t-τ)𝐅(n,τ) transitions the system in the interval [τ,t] to a final probability vector, during which no total change to n is made. This contributes to the probability vector 𝐏(n,t); the final step integrates over all possible first-passage times:


    𝐏(n,t)     = ∫^t_0 dτ𝐓(0,t-τ)𝐅(n,τ).


From here it is easier to work in Laplace space. The RHS of Eq.(<ref>) is a convolution, and using Eq.(<ref>) one obtains 


    𝐓̃(n|z)𝐏(0)     = 𝐓̃(0,z)𝐅̃(n,z).


In Laplace space, the total first-passage time probability is 


    F̃(n|z)     = (𝐈,𝐓̃(0,z)^-1𝐓̃(n|z)𝐏(0)) .


As Ptaszynski <cit.> notes, in order to connect the FPTD to the WTD, the initial state 𝐏(0) needs to be the normalised system state just after a forward jump to the drain has occurred:


    𝐏(0)     = 𝐉_F𝐏̅/(𝐈,𝐉_F𝐏̅).

Combining this choice of 𝐏(0) with Eq.(<ref>), we obtain 


    F̃(n|z)     = (𝐈,𝐓̃(0,z)^-1𝐓̃(n|z)𝐉_F𝐏̅)/(𝐈,𝐉_F𝐏̅),

where F(n|τ), defined as the inverse Laplace transform of Eq.(<ref>) ℒ^-1{F̃(n|z)}, is the probability that, given an initial tunneling to the drain occurred at some arbitrary time in the stationary state, the jump number first reached n after a time interval τ. For unidirectional transport, F(1|τ) exactly reproduces the WTD defined in Eq.(<ref>).

The FPTD in time-space is obtained from the numerical inverse Laplace transform of Eq.(<ref>):


    F(n|τ)     = 1/2π ilim_R→∞∫_c-iR^c+iRdz   e^zτ(𝐈,𝐓̃(0,z)^-1𝐓̃(n|z)𝐉_F𝐏̅)/(𝐈,𝐉_F𝐏̅).


We focus, however, on the cumulants of the FPTD, which are again calculated from F̃(n|z):


    ⟨⟨τ^k_n⟩⟩^*    = (-1)^klim_z→ 0^+[d^k/dz^klnF̃(n|z)],

where the limit z→ 0^+ is required as 𝐋(χ) is singular for χ = {0,2π} <cit.>. We have introduced the ^* notation for FPTD cumulants to distinguish them from WTD cumulants, and used Ptaszynski's notation of ⟨⟨τ^k_n⟩⟩^* referring to the k^th cumulant of F(n|τ).

Ptaszynski has recently shown that when the transport is renewal the FPTD cumulants are linearly related:

    ⟨⟨τ^k_n⟩⟩^*    = n⟨⟨τ^k_1⟩⟩^*,

as well as deriving, similarly to the WTD, renewal relations between the FCS and the first-passage time cumulants <cit.>. The average current is then 


    ⟨ I ⟩    = n/⟨τ_n⟩^*
        = 1/⟨τ_1⟩^*.

Unlike the WTD, where the equality in Eq.(<ref>) is satisfied regardless of renewal behaviour, the equality in Eq.(<ref>) is only true for renewal processes; the logic applied in Eq.(<ref>) is not available to first-passage times since, for non-renewal processes, F(2|τ) ≠ F(1|τ;1|τ').

One can also define a randomness parameter R^*_n for each distribution F(n|τ):


    R^*_n    = n⟨⟨τ_n^2⟩⟩^*/(⟨τ_n⟩^*)^2.

For renewal systems they all reduce to R^*_1, which is then equal to the Fano factor. Unlike R, we expect that the renewal relationships will hold in bidirectional transport. The Pearson correlation coefficient between subsequent first-passage times can be calculated using F(2|τ) <cit.>:


    p^*    = 2R^*_2/R^*_1 - 1.


Throughout the majority of the paper we will not use the subscript n for the cumulants, with the implication being that we are referring to cumulants of F(1|τ) unless explicitly stated.





§ MODEL SYSTEMS
 

In this section we introduce the three transport scenarios under consideration and define the rate equations describing their dynamics. Quantum rates are calculated using the T-matrix method, which is a suitable replacement for more rigorous approaches when the transport is Markovian; that is, when T_α,V_SD≫γ, where γ is the broadening of the molecular energy levels, T_α is the temperature of electrode α, and V_SD is the source-drain voltage bias <cit.>. In order to use rate equations as opposed to comprehensive master equations <cit.> we also operate under the secular approximation; coherence between molecular basis states are ignored. All scenarios describe the transport of electrons through a molecule weakly coupled to two macroscopic metal electrodes. The total Hamiltonian is 


    H = H_S+H_D+H_M+H_T,

where H_M is the Hamiltonian of the nanoscale quantum system; H_S and H_D are the source and drain Hamiltonians

    H_S+H_D = ∑_α = S,D∑_k,σε^_α,k,σa^†_α,k,σa^_α,k,σ;

and H_T is the system-electrode interaction

    H_T = ∑_α = S,D∑_k, σt^_α,k,σ(a^†_α,k,σa^_σ+a_σ^†a^_α,k,σ).

The operator a_α,k,σ^† creates an electron in electrode α in the single-particle state k with spin σ = {↑,↓} and free energy ε_α,k,σ, while a_α,k,σ annihilates an electron with spin σ from the same state in electrode α, and t_α,k,σ are the tunneling amplitudes between the quantum system and the electrodes.

The molecule is described by the Anderson-Holstein model; an energy level, which can be spin split, is coupled to a vibrational mode such that electrons occupying the molecule experience electron-electron and electron-phonon interactions.


    H_M    = ∑_σε_σ a_σ^†a_σ+U n_↑n_↓ + λ(b^† + b)n+ω b^†b.


The operators a^†_σ and a_σ create and annihilate an electron in the single-particle state σ with free energy ε_σ, while b^† and b create and annihilate phonons with frequency ω. The particle number operators are n_σ = a^†_σa_σ and n = ∑_σ a^†_σa_σ. The λ(b^† + b)n term describes the electron-phonon interaction, where λ is the coupling strength. Likewise, U is the Coulomb repulsion strength. 

The transition rate from the electrode-system many-body state i to many-body state f is 


    Γ_fi = 2π|⟨f|T|i⟩|^2δ(E_i-E_f),
 where the T-matrix is 

    T     = H_T+H_T1/E_i-H_0+iηH_T
        + H_T1/E_i-H_0+iηH_T1/E_i-H_0+iηH_T+... .


The leading order approximation is Fermi's golden rule, and describes sequential tunneling:

    Γ_fi = 2π|⟨f|H_T|i⟩|^2δ(E_i-E_f).

Fermi's generalised golden rule <cit.> details the higher order terms; for example, cotunneling events are described by going to fourth order in H_T. 



 §.§ Anderson impurity


For a rigid molecule with no electron-phonon coupling λ = 0, Eq.(<ref>) reduces to the Anderson model:


    H_M    = ∑_σε_σ a^†_σa_σ + U n_↑n_↓̅.


The molecular orbital is then in one of four system states; it can be empty |0⟩, occupied by a single electron of either spin |σ⟩, or occupied by two electrons of opposite spin |2⟩, with probability vector 𝐏(χ,t) = [P_0(χ,t),P_↑(χ,t),P_↓(χ,t),P_2(χ,t)]^T. Including cotunneling and sequential tunneling events, the Liouvillian minus the forward and backward jump operators is 


    𝐋_0    =  [[              -(Γ_↑0+Γ_↓0 +Γ_00^(2))                              Γ_0↑^S                              Γ_0↓^S                                   0;                                    ;                             Γ_↑ 0^S    -(Γ_0↑+Γ_2↑ + Γ_↑↑^(2)+Γ_↓↑^(2))                     Γ_↑↓^SS+Γ_↑↓^DD                              Γ_↑2^S;                                    ;                              Γ_↓0^S                     Γ_↓↑^SS+Γ_↓↑^DD -(Γ_0↓+Γ_2↓ + Γ_↓↓^(2) +  Γ_↑↓^(2))                              Γ_↓2^S;                                    ;                                   0                              Γ_2↑^S                              Γ_2↓^S              -(Γ_↑2+Γ_↓2 +Γ_22^(2)) ]],


with the forward and backward jump operators defined as 



    𝐉_F    = [[ Γ_00^SD  Γ_0↑^D  Γ_0↓^D       0;                ;       0 Γ_↑↑^SD Γ_↑↓^SD  Γ_↑2^D;                ;       0 Γ_↓↑^SD Γ_↓↓^SD  Γ_↓2^D;                ;       0       0       0 Γ_22^SD ]]      and    𝐉_B = [[ Γ_00^DS       0       0       0;                ;  Γ_↑0^D Γ_↑↑^DS Γ_↑↓^DS       0;                ;  Γ_↓0^D Γ_↓↑^DS Γ_↓↓^DS       0;                ;       0  Γ_2↑^D  Γ_2↓^D Γ_22^DS ]].


The total sequential rates are 

    Γ_kl    = Γ_kl^S + Γ_kl^D;

the total elastic cotunneling rates are 

    Γ_ll^(2)    = Γ_ll^SD + Γ_ll^DS;

and the total inelastic cotunneling rates are 

    Γ_kl^(2)    = Γ_kl^SD + Γ_kl^DS + Γ_kl^SS + Γ_lk^DD.


The sequential rates are


    Γ_σ0^α    = γ^α  n_F(ε_σ-μ_α), 
    Γ_0σ^α    = γ^α (1-n_F(ε_σ-μ_α)), 
    Γ_σ2^α    = γ^α (1-n_F(ε_σ+U-μ_α)),    and
    Γ_2σ^α    = γ^α  n_F(ε_σ+U-μ_α),

where γ^α=2π|t_α|^2ρ(ε_α); ρ(ε_α) is the density of states for the source and drain electrodes, which is assumed to be constant; and n_F(ε-μ_α) is the Fermi-Dirac distribution:

    n_F(ε-μ_α)      = 1/1+e^(ε-μ_α)/T_α.


During cotunneling, the electron undergoes two tunneling events in the same quantum process through an intermediate “virtual” energy level in the molecule. Elastic cotunneling leaves the molecule in the same energetic state as before the process, while inelastic cotunneling leaves the molecular orbital occupied by the same number of electrons but in an excited or de-excited state. The notation Γ_kl^αβ implies that an electron cotunneled from electrode α to electrode β, changing the system from state l to state then k; if the process is elastic l=k. Elastic cotunneling processes do not appear in Eq.(<ref>) as they do not affect the molecular probabilites 𝐏(t), while elastic cotunneling events involving the same electrode α=β are not included even in the n-resolved master equation as they neither change the jump number or the molecular state.

The cotunneling rates are calculated by going to fourth order in H_T in Eq.(<ref>) and Eq.(<ref>). To avoid divergences in the T-matrix calculated cotunneling rates, and overcounting of sequential processes, one must use the regularisation procedure noted first by Averin <cit.> and detailed separately by Koch et al. <cit.> and Turek and Matveev <cit.>. The details and specific form of the lengthy cotunneling rates for an Anderson impurity are in Appendix <ref>.



 §.§ Holstein model


In the limit U→∞, the molecule cannot be doubly occupied and behaves as a single level coupled to a vibrational mode: the Holstein model, with the Hamiltonian


    H_M    = ε_0 a^†a + λ(b^† + b)a^†a+ω b^†b.


After the Lang-Firsov transformation <cit.> the molecular Hamiltonian reduces to 

    H_M    = εã^†ã + ωb̃^†b̃,

with ε = ε_0 - λ^2/ω and new fermionic ã^†ã and bosonic b̃^†b̃ operators obtained from the original operators via the canonical transformation. The molecular eigenstate is now described by two quantum numbers |nq⟩, where n={0,1} is the electron occupation number and q = {0,1,...,+∞} is the vibrational quanta occupation number. The associated eigenenergy is E_nq = ε n + ω q. The Fourier transformed vector is then

    𝐏(χ,t)     = [[ P_00(χ,t); P_10(χ,t); P_01(χ,t); P_11(χ,t);         .;         .;         .; P_0N(χ,t); P_1N(χ,t) ]],

which is described by the master equation <cit.>





    Ṗ_0q(χ,t)     = ∑_q'(Γ^S_0q;1q' + Γ^D_0q;1q'e^iχ)P_1q'(χ,t) - ∑_α q'Γ^α_1q';0qP_0q(χ,t), 
    Ṗ_1q(χ,t)     = ∑_q'(Γ^S_1q;0q' + Γ^D_1q;0q'e^-iχ)P_0q'(χ,t) - ∑_α q'Γ^α_0q';1qP_1q(χ,t) .


Evidently, there can be any number of interacting vibrations and 𝐏(χ,t) is a vector of infinite length. In practice, however, we only include a finite number N of vibrations in the transport calculations, so that 𝐏(χ,t) and 𝐈 are vectors of length 2(N+1). Likewise, the Liouvillian and two jump operators are 2(N+1)×2(N+1) matrices. 𝐉_F and 𝐉_B are explicitly written in Appendix <ref>.

The rates are 

    Γ^α_1q';0q    = γ|X_q'q|^2n_F(ε + ω(q'-q) - μ_α)  and 
    Γ^α_0q';1q    = γ|X_q'q|^2[1-n_F(ε + ω(q'-q) - μ_α)],


where the Franck-Condon matrix elements are 

    X_qq'    = ⟨ q | e^-λ(b^†-b)|q'⟩.


We consider two distinct vibrational scenarios. So far we have described transport when the phonons are unequilibrated, but when the phonons are in equilibrium with an external bath at temperature T_V, the rate equation reduces to that of a single resonant level, with Liouvillian


    𝐋_0    = [[  -T_10 T_01^S; T^S_10  -T_01 ]],
 and jump operators



    𝐉_F    = [[      0 T_01^D;      0      0 ]],    𝐉_B    = [[      0      0; T_10^D      0 ]],


where the probabilities are defined with the ansatz 

    P_nq(χ,t)     = P_n(χ,t)e^-qω/T_V/1-e^-ω/T_V;

the total transition rates are T_kl = ∑_α T_kl^α; and the individual transition rates are 

    T^α_kl = ∑_qq'Γ^α_kq;lq'e^-qω/T_V/1-e^-ω/T_V.






§ RESULTS
 

In this section we compare the full counting statistics, waiting time distribution, and first-passage time distribution for the three transport scenarios.



 §.§ Sequential tunneling and cotunneling through an Anderson impurity


In a previous paper, we examined, via waiting times, the effect cotunneling has on strong electron-electron interactions in an Anderson impurity <cit.>. In order to include elastic cotunneling events, it is necessary to use the n-resolved master equation to calculate the waiting times. Unfortunately, when bidirectional transitions are included the definition of the idle-time probability in Eq.(<ref>) is no longer satisfied <cit.>. One may of course treat the transport as unidirectional and calculate the waiting times in either direction; Fig.(<ref>) shows that doing this, in a regime where backtunneling processes are present, does not reproduce the correct cotunneling current from 1/⟨τ⟩_F,co - 1/⟨τ⟩_B,co. The first-passage time prediction 1/⟨τ⟩_co^*, on the other hand, reproduces the cotunneling current; the black and red dashed lines in Fig.(<ref>) coincide. Over the voltage range 2-3meV there is a small numeric difference between ⟨ I ⟩_co and 1/⟨τ⟩_co^*, since the transport is non-renewal, but it is not visible to the naked eye.  In comparison, the sequential current is exactly reproduced by 1/⟨τ⟩_F,seq - 1/⟨τ⟩_B,seq, since sequential tunneling is described completely by a standard master equation, and 1/⟨τ⟩_seq^*, since the sequential transport is always renewal. One can see this in Fig.(<ref>), as all solid lines coincide.



These difficulties previously kept us to the large bias voltage regime. It is well known in the literature, however, that in a certain voltage range inelastic cotunneling processes cause telegraphic switching in an Anderson impurity, which is now accessible to fluctuating time statistics via the FPTD. The Fano factor in Fig.(<ref>) is a reproduction of Fig.(3) from Ref.[Thielmann2005a], with the first-passage and waiting time randomness parameters included as well. The molecule is in the Coulomb blockade regime U≫ V_SD, so that transport through one σ level blocks transport through the other σ̅ level. At small voltages, for sequential tunneling, Belzig <cit.> has shown that the molecule goes through cycles of transferring n-1 electrons quickly through the ↑ level until it gets “stuck” in the ↓ level. The cycle completes when the ↓ empties via thermal broadening of the drain occupation function. The noise is thus characteristic of a sum of Poissonian processes, each transferring n electrons, and weighted by (1/2)^n. In contrast, when cotunneling processes are included, super-Poissonian noise arises from a different mechanism. Note that from here we have not included cumulants from the forward cotunneling WTD, such as R and p, since in this voltage regime the transport is bidirectional. 

At low bias V_SD/2 < ε_↑ the only transport mechanism is elastic cotunneling through the empty ↑ level or the occupied ↓ level; the noise is Poissonian and F=1. When ε_↓ > V_SD/2 ≥ε_↑ the ↑ level may be occupied via sequential tunneling, and there are now two main transport channels: elastic cotunneling through the ↓ level and sequential tunneling through ↑ level. Telegraphic switching between these two processes occurs when inelastic cotunneling induces a spin flip σ→σ̅ or the molecule relaxes from ↑→↓ via sequential tunneling. Super-Poissonian noise is a symptom of the switching between two Poissonian processes, which is the main observation of Fig.(3) from Ref.[Thielmann2005a] and is evident in Fig.(<ref>). In this voltage regime the transport is evidently non-renewal as R^* does not match the Fano factor; there should thus be correlations between subsequent first-passage times. 

Positive correlations in the waiting times are known to occur for systems undergoing telegraphic switchin: for example, in systems with spin-polarised electrodes <cit.>. The current switches between two transport channels, each with a characteristic waiting time. When in the faster channel, waiting times between successive tunnelings are more likely to be shorter, and vice versa in the slower channel. Consequently, waiting times between successive tunnelings to the drain are positively correlated; one would expect that there would be similar positive correlations between the first-passage times. In Fig.(<ref>), however, we report intriguing negative correlations between successive first-passage times. That is, if the jump number reaches +1 after a short time-interval, it is more likely that the time until the jump number next reaches +1 will be long. The values are small enough, however, that the first-passage times are still not considered correlated. A possible explanation is that the time it takes for the transport to switch between the two mechanisms is comparable to the average first-passage times for the two channels, so that the system does not spend long enough in either channel for subsequent first-passage times to become positively correlated.







 §.§ Vibrations with a finite bias voltage


One of the authors (Kosov) has previously reported positive correlations in successive waiting times for transport through the Holstein model <cit.>, which we now compare to first-passage times.

We consider transport through the Holstein model with fixed temperature and over a range of V_SD. In such a scenario, Koch and von Oppen have previously explained the different behaviour of equilibrated and unequilibrated current at low voltages <cit.>, which originates from X_qq' being suppressed on diagonal and close-to-diagonal elements under a strong electron-phonon coupling λ≫ 1. At small T and eV, when the electrons are not energetic enough to reach highly excited phonon states, the current is thus suppressed, which is evident more for equilibrated phonons than out-of-equilibrium phonons, as multiple unequilibrated transitions can produce the required excited phonon states. At higher voltages, the current also displays the characteristic steplike structure, with steps occurring when the voltage is a multiple of 2, since ω = 1.

The total electronic current is correctly predicted from the WTD in Eq.(<ref>) for both equilibrated and unequilibrated phonons; the blue and black dashed and solid lines in Fig.(<ref>) exactly overlap. In comparison, 1/⟨τ⟩^* correctly predicts the current over all voltages for equilibrated but not unequilibrated phonons, implying that at the voltages where they differ the transport is non-renewal. One can see that the red dashed line in Fig.(<ref>) coincides with both the blue and black dashed lines, but that there is a difference between the red solid line and the blue and black solid lines. Unlike the current plots for the other transport scenarios we consider, the discrepancy between the unequilibrated 1/⟨τ⟩^* and ⟨ I ⟩ due to renewal processes is large enough to be visually recognisable. As Fig.(<ref>) shows, however, the largest correlations occur at voltages where the current, and thus the difference between 1/⟨τ⟩^* and ⟨ I ⟩, is negligible; consequently, we turn to the randomness parameters to analyse the non-renewal behaviour.

As found in Ref.[Kosov2017b], when phonons are unequilibrated, there exist strong correlations between successive waiting times, due to an elastic “shortcut” channel that opens in a small voltage range. In Fig.(<ref>) we report similar correlations in the first-passage times over the same voltage range, albeit with a smaller peak. This is expected; due to backtunneling events the first-passage time can never be smaller than the waiting time τ^*≥τ, and so correlations between first-passage times must be smaller. The Fano factor and randomness parameters in Fig.(<ref>) support this non-renewal behaviour. At low voltages R fails to capture the divergence of the Fano factor and R^*. Between the 1→2eV voltage range F and R^* differ, correctly predicting the non-renewal behaviour in Fig.(<ref>). In comparison, R actually crosses the F at ∼ 1.9eV, implying renewal behaviour, while at the same voltage the waiting time Pearson coefficient is still a non-zero p ≈ 0.05. The correlations disappear when the phonons are forced to relax to equilibrium immediately, as the excited phonon state providing the elastic shortcut channel is no longer occupied. 







 §.§ Vibrations with a finite temperature gradient
 

At finite temperatures, and when the electrode Fermi energies are approximately in resonance with ε + ω q, the bidirectional transport is dominated by thermal noise. A recent experimental paper has shown, however, that electronic noise distinct from the thermal noise can be generated from a temperature gradient Δ T in the absence of a bias voltage, which they label delta-T noise <cit.>. We finish our investigation by exploring the effect a temperature gradient has on the first-passage times for sequential tunneling through the Holstein model.

Since the transport in this scenario is described by a standard master equation, the exact current matches the prediction from the WTD in Eq.(<ref>), which is shown in Fig.(<ref>); all black and blue lines coincide. The FPTD current prediction from Eq.(<ref>), on the other hand, differs from ⟨ I ⟩ at low Δ T, implying non-renewal behaviour in this regime; otherwise, the red FPTD lines coincide with all blue and black lines. Interestingly, the equilibrium current is larger than the non-equilibrium current, which is understood by examining X_qq' and the phonon occupation probabilities. For λ = 1 and ω = 1 the Franck-Condon factor is largest for elastic q=0 transitions. These transitions are more accessible to phonons in equilibrium, as the q=0 state has a larger relative occupation probability for equilibrium phonons than non-equilibrium phonons: 


    P^eq_q=0/P^eq_q=1    > P^ne_0;q=0+P^ne_1;q=0/P^ne_0;q=1+P^ne_1;q=1.

We note that for large Δ T the total current is the same order as that induced by a voltage bias in Fig.(<ref>).

In Fig.(<ref>) F diverges as Δ T → 0, which is expected since thermal processes dominate in this regime. At larger Δ T the noise no longer diverges but is still super-Poissonian, and distinct from the shot noise generated by a voltage bias. We attribute this to the relatively large temperatures; even the smallest drain temperature we consider is large enough to cause backtunnelings from the drain, hence the super-Poissonian partition noise normally associated with a voltage bias.

It is evident that R^*≠ F at small Δ T, implying that the transport is non-renewal in this temperature range. Despite this non-renewal behaviour, we have chosen not to display the Pearson correlations: they are non-zero but small. At larger Δ T when the renewal assumption is satisfied, R^* successfully predicts F, while R clearly does not as it only includes unidirectional transitions.



§ CONCLUSION
 

In this paper we have compared the waiting time statistics to the first-passage time statistics for a molecule with strong electron-electron or electron-phonon interactions. All methods rely on a form of the master equation, in which the on-site interactions are treated exactly and the coupling between the molecule and electrode is treated perturbatively under the Born-Markov approximation. In bidirectional transport scenarios we have demonstrated the effectiveness of using first-passage times as opposed to waiting times. In particular, we compared the total current to the current computed from the WTD and FPTD, the Fano factor to the WTD and FPTD randomness parameters, and the correlation between subsequent waiting and first-passage times. We first examine sequential tunneling and cotunneling through an Anderson impurity. When a full transport description requires use of the n-resolved master equation, the WTD can be calculated in the forward tunneling regime only. Such difficulties arise when including elastic cotunneling, but do not apply to first-passage times; we demonstrate that there are unexpected, minor negative correlations between the first-passage times even though the transport displays telegraphic switching. Next, we analysed the Holstein model at finite bias voltage, in which correlations between successive first-passage times are smaller than correlations between successive waiting times when backtunneling events are present. Finally, we considered transport through the Holstein model at zero bias voltage and a finite temperature gradient. The temperature gradient induces a current, but also significant backtunneling events. Hence, the forward waiting time τ_F was incapable of accurately predicting the Fano factor, as opposed to the first-passage time τ^*.





§ COTUNNELING RATES
 

Cotunneling rates calculated from the T-matrix approach are regularised according to the procedure outlined by Koch et al <cit.>. Calculating cotunneling rates for transport through an Anderson impurity, in particular, has been thoroughly covered by in Ref.[Kaasbjerg2015,Rudge2018], to which we direct the reader for a comprehensive overview.

The elastic cotunneling rates are 



    Γ_nn^αβ = γ^2/2πlim_η→0^+∫ dε| 1/ε-E_1+iη±1/ε-E_2+iη| n_F(ε-μ_α)[1-n_F(ε-μ_β)],


where E_1 and E_2 are derived from the cotunneling pathways involved in the rate, and the ± is only positive for elastic tunneling through an initially empty or initially doubly occupied system. Similarly, inelastic cotunneling rates have the general form 


    Γ_σ̅σ^αβ = γ^2/2πlim_η→0^+∫ dε| 1/ε-ε_σ̅-U+iη - 1/ε-ε_σ̅-iη| n_F(ε-μ_α)[1-n_F(ε-μ_β+ε_σ-ε_σ̅)],


where the broadening of the intermediate virtual state iη in Eq.(<ref>) and Eq.(<ref>)was first noted by Averin. 

The rates are then split into a sum via:



    | 1/ε-E_1+iη±1/ε-E_2-iη|     = 1/(ε-E_1)^2+η^2+1/(ε-E_2)^2+η^2±2{1/ε-E_1+iη·1/ε-E_2-iη}.



All that remains is to compute the two types of integrals:


    lim_η→ 0∫ dεn_F(ε - μ_α)[1-n_F(ε-μ_β)]/(ε - E_1-iη)·(ε - E_2+iη)    = 
    n_B(μ_β-μ_α)/E_2-E_1{ψ^(1)(E^α_-,2) - ψ^(1)(E^α_+,1) 
                        - ψ^(1)(E^β_-,2) + ψ^(1)(E^β_+,1)}, 
    lim_η→ 0∫ dεn_F(ε - μ_α)[1-n_F(ε-μ_β)]/(ε - E_k)^2+η^2 - 𝒪(η^-1)     = 
    n_B(μ_β-μ_α)/2π T{ψ^(1)( E^β_+,k) - ψ^(1)( E^α_+,k)}.



Here, E^α_±,k = 1/2 ±  iμ_α - E_k/2π T; ψ(x) and ψ^(1)(x) are the digamma and trigamma functions, respectively; and the Bose-Einstein distribution is 


    n_B(μ_β-μ_α)     = 1/e^(μ_β-μ_α)/T-1.
 



§ JUMP OPERATORS FOR THE HOLSTEIN MODEL
 

The jump operators defined from the master equation in Eq.(<ref>) and Eq.(<ref>) are 



    𝐉_F    = [[             0     Γ^D_00;10             0     Γ^D_00;11             0             ⋯             ⋯             ⋯     Γ^D_00;1N;             0             0             0             0             0                                                       0;             0     Γ^D_01;10             0     Γ^D_01;11             0                                               Γ^D_01;1N;             ⋮                                         ⋱                                                                     ⋮;             ⋮                                                       ⋱                                                       ⋮;             ⋮                                                                     ⋱                                         ⋮;             0                                                                     0             0             0             0;             0                                                                     0 Γ^D_0N;1(N-1)             0     Γ^D_0N;1N;             0             ⋯             ⋯             ⋯             ⋯             0             0             0             0;               ]],  and 
    𝐉_B    = [[             0             0             0             0             ⋯             ⋯             ⋯             0             0;     Γ^D_10;00             0     Γ^D_10;01             0                                               Γ^D_10;0N             0;             0             0             0             0                                                       0             0;     Γ^D_11;00             0     Γ^D_11;01             0                                               Γ^D_11;0N             0;             ⋮                                                       ⋱                                                       ⋮;             ⋮                                                                     ⋱                                         ⋮;             ⋮                                                                                   ⋱                           ⋮;             0                                                                     0             0             0             0;     Γ^D_1N;00             ⋯             ⋯             ⋯             ⋯ Γ^D_1N;0(N-1)             0     Γ^D_1N;0N             0;               ]],






