Non-renewal statistics for electron transport in a molecular junction with electron-vibration interaction

TL;DR

This paper investigates how electron-vibration interactions in a molecular junction lead to non-renewal statistics in electron tunnelling times, revealing correlations influenced by vibrational dynamics and coupling strength.

Contribution

It introduces a master equation approach to analyze non-renewal statistics in electron transport, highlighting the role of vibrational dynamics and strong electron-vibration coupling.

Findings

Non-renewal statistics emerge when vibrational dynamics are treated exactly.

Strong electron-vibration coupling is necessary for temporal correlations.

Rare tunnelling events in Franck-Condon blockade are positively correlated.

Abstract

Quantum transport of electrons through a molecule is a series of individual electron tunnelling events separated by stochastic waiting time intervals. We study the emergence of temporal correlations between successive waiting times for the electron transport in a vibrating molecular junction. Using master equation approach, we compute joint probability distribution for waiting times of two successive tunnelling events. We show that the probability distribution is completely reset after each tunnelling event if molecular vibrations are thermally equilibrated. If we treat vibrational dynamics exactly without imposing the equilibration constraint, the statistics of electron tunnelling events become non-renewal. Non-renewal statistics between two waiting times $\tau_1$ and $\tau_2$ means that the density matrix of the molecule is not fully renewed after time $\tau_1$ and the probability of observing waiting time $\tau_2$ for the second electron transfer depends on the previous electron waiting time $\tau_1$. The strong electron-vibration coupling is required for the emergence of the non-renewal statistics. We show that in Franck-Condon blockade regime the extremely rare tunnelling events become positively correlated.

Equations

$$H=H_{\text{molecule}}+H_{\text{electrodes}}+H_T.$$

$$H_{\text{molecule}}={ϵ}_0{a}^{†}a+\lambda \omega (b^{†}+b)a^{†}a+\omega b^{†}b,$$

$$H_{\text{electrodes}}=\sum _{k,\alpha =S,D}{ϵ}_{k\alpha }{a}_{k\alpha }^{†}{a}_{k\alpha },$$

$$H_T=\sum _{k,\alpha =S,D}{t}_{\alpha }(a_{k\alpha }^{†}a+a^{†}{a}_{k\alpha }),$$

$${\dot{P}}_{0q}(t)$$

$${\dot{P}}_{1q}(t)$$

$${\Gamma }_{0q^{\prime },1q}^{\alpha }={\gamma }^{\alpha }|X_{q^{\prime }q}{|}^2\left(1-f_{\alpha }[ϵ-\omega (q^{\prime }-q)]\right).$$

$${\Gamma }_{1q^{\prime },0q}^{\alpha }={\gamma }^{\alpha }|X_{q^{\prime }q}{|}^2{f}_{\alpha }[ϵ+\omega (q^{\prime }-q)]$$

$$f_{\alpha }(E)=\frac{1}{1+e^{(E-{\mu }_{\alpha })/T}},$$

$$X_{q{q}^{\prime }}=⟨q|e^{-\lambda (b^{†}-b)}|q^{\prime }⟩,$$

$${\gamma }^{\alpha }=2\pi t_{\alpha }^2{\rho }_{\alpha }$$

$$\mathbf{P}(t)=\left[\begin{array}{c}{P}_{00}(t)\\ P_{10}(t)\\ P_{01}(t)\\ P_{11}(t)\\ ⋮\\ P_{0N}(t)\\ P_{1N}(t)\end{array}\right],$$

$$\mathbf{I}=\left[\begin{array}{c}1\\ 1\\ 1\\ 1\\ ⋮\\ 1\\ 1\end{array}\right].$$

$$(\mathbf{I},\mathbf{P}(t))=\sum _{q=0}^N{P}_{0q}(t)+P_{1q}(t)=1.$$

$$\dot{\mathbf{P}}(t)=\mathcal{L}\mathbf{P}(t),$$

$$(J\mathbf{P}(t){)}_{mq}={\delta }_{m0}\sum _{q^{\prime }}{\Gamma }_{0q,1q^{\prime }}^D{P}_{1q^{\prime }}(t).$$

$$\mathcal{L}\mathbf{P}=0.$$

$$w_2({\tau }_2,{\tau }_1)=(\mathbf{I},J{e}^{(\mathcal{L}-J){\tau }_2}J{e}^{(\mathcal{L}-J){\tau }_1}J\mathbf{P}).$$

$$w_1(\tau )=(\mathbf{I},J{e}^{(\mathcal{L}-J)\tau }J\mathbf{P}).$$

$${\int }_0^{\infty }d{\tau }_1{\int }_0^{\infty }d{\tau }_2{w}_2({\tau }_2,{\tau }_1)$$

$$=(\mathbf{I},J(\mathcal{L}-J{)}^{-1}J(\mathcal{L}-J{)}^{-1}J\mathbf{P})$$

$$=(\mathbf{I},(J-\mathcal{L}+\mathcal{L})(\mathcal{L}-J{)}^{-1}J(\mathcal{L}-J{)}^{-1}J\mathbf{P})$$

$$=-(\mathbf{I},J(\mathcal{L}-J{)}^{-1}J\mathbf{P})$$

$$=-(\mathbf{I},(J-\mathcal{L}+\mathcal{L})(\mathcal{L}-J{)}^{-1}J\mathbf{P})$$

$$=(\mathbf{I},J\mathbf{P}).$$

$$w_2({\tau }_2,{\tau }_1)=\frac{(\mathbf{I},J{e}^{(\mathcal{L}-J){\tau }_2}J{e}^{(\mathcal{L}-J){\tau }_1}J\mathbf{P})}{(\mathbf{I},J\mathbf{P})}.$$

$$w_1(\tau )=\frac{(\mathbf{I},J{e}^{(\mathcal{L}-J)\tau }J\mathbf{P})}{(\mathbf{I},J\mathbf{P})}.$$

$${\int }_0^{\infty }d{\tau }_2{w}_2({\tau }_2,{\tau }_1)$$

$$=-\frac{(\mathbf{I},J(\mathcal{L}-J{)}^{-1}J{e}^{(\mathcal{L}-J){\tau }_1}J\mathbf{P})}{(\mathbf{I},J\mathbf{P})}$$

$$=-\frac{(\mathbf{I},(J-\mathcal{L}+\mathcal{L})(\mathcal{L}-J{)}^{-1}J{e}^{(\mathcal{L}-J){\tau }_1}J\mathbf{P})}{(\mathbf{I},J\mathbf{P})}$$

$$=\frac{(\mathbf{I},J{e}^{(\mathcal{L}-J){\tau }_1}J\mathbf{P})}{(\mathbf{I},J\mathbf{P})}=w_1({\tau }_1).$$

$${\int }_0^{\infty }d{\tau }_1{w}_2({\tau }_2,{\tau }_1)$$

$$=-\frac{(\mathbf{I},J{e}^{(\mathcal{L}-J){\tau }_2}J(\mathcal{L}-J{)}^{-1}J\mathbf{P})}{(\mathbf{I},J\mathbf{P})}$$

$$=-\frac{(\mathbf{I},J{e}^{(\mathcal{L}-J){\tau }_2}J(\mathcal{L}-J{)}^{-1}(J-\mathcal{L}+\mathcal{L})\mathbf{P})}{(\mathbf{I},J\mathbf{P})}$$

$$=\frac{(\mathbf{I},J{e}^{(\mathcal{L}-J){\tau }_2}J\mathbf{P})}{(\mathbf{I},J\mathbf{P})}=w_1({\tau }_2).$$

$$K(x_1,x_2)={\int }_0^{\infty }d{\tau }_1{\int }_0^{\infty }d{\tau }_2{e}^{i{x}_1{\tau }_1}{e}^{i{x}_2{\tau }_2}{w}_2({\tau }_1,{\tau }_2).$$

Full text

I Introduction

On the microscopic quantum mechanical level, electron current consists of a sequence of single electron tunneling events separated by random waiting time intervals. Nazarov and Blanter (2009) Statistics of these waiting time intervals reveals the wealth of interesting information about details of quantum transport. Thomas and Flindt (2013); Sothmann (2014); Thomas and Flindt (2014); Potanina and Flindt (2017); Tang et al. (2014); Seoane Souto et al. (2015); Goswami and Harbola (2015); Rudge and Kosov (2016a, b); Kosov (2017) The statistical properties of the waiting times are usually studied using waiting time distribution (WTD), which is a conditional probability distribution that we observe the electron transfer in the detector electrode (drain or source) at time $t+\tau$ given that an electron was detected in the same electrode at time $t$.Brandes (2008) WTD is a complementary to very popular full counting statistics in quantum transport and it has recently gained a significant popularity in the study of nanoscale and mesoscale systems.Albert et al. (2012); Thomas and Flindt (2013); Sothmann (2014); Thomas and Flindt (2014); Dasenbrook et al. (2015); Tang et al. (2014); Seoane Souto et al. (2015); Goswami and Harbola (2015); Rudge and Kosov (2016a, b); Kosov (2017)

The question which we discuss in this paper is the following. When an electron transfers through a molecular junction, is there exists a correlation between waiting times for successive electron tunneling or they are statistically independent? Is it possible to have, for example, a situation, when the second tunneling electron senses the waiting time of the previous electron and changes its own waiting time accordingly? Intuitively, we expect that such kind of statistical temporal correlations can emerge in molecular junctions with strong electron-vibrational coupling – the excitation of a particular vibrational state depends on the waiting time of the electron, and the waiting time for the second electron feels the previous one through this vibration. Strong coupling between electronic and nuclear dynamics distinguishes molecular junctions from other nanoscale quantum transport systems. The interplay between nuclear and electronic dynamics has already led to the discovery of distinctly molecular junction phenomena as Franck-Condon blockade, Koch and von Oppen (2005); Koch et al. (2006); Lau et al. (2016) negative differential resistance, Härtle and Thoss (2011); Kuznetsov (2007); Galperin et al. (2005); Zazunov et al. (2006); Dzhioev and Kosov (2012) nonequilibrium chemical reactions,Dzhioev and Kosov (2011); Thomas et al. (2012); Dzhioev et al. (2013) cooling of nuclear motion by electric current.Galperin et al. (2009); Ioffe et al. (2008); Härtle and Thoss (2011)

The basis for our work is the extension of ideas of WTD beyond studying statistics of single waiting time to the domain of multiple waiting times joint probability distribution.Ptaszyński (2017) Renewal theory assumes that successive waiting times between transport events are statistically independent equally distributed random variables. In this case the joint probability density of two successive waiting times $w_{2}(\tau_{2},\tau_{1})$ can be factorised into a product of two single-time distributions $w_{1}(\tau_{2})w_{1}(\tau_{1})$, that means that the distribution is totally ”renewed” after waiting time $\tau_{1}$. The non-renewal statistics means the existence of temporal correlations between the subsequent tunneling events $w_{2}(\tau_{2},\tau_{1})\neq w_{1}(\tau_{2})w_{1}(\tau_{1})$. The study of non-renewal statistics is a very interesting topic on its own. Although the questions of temporal correlations between electron tunneling events are only started to appear in quantum transport,Ptaszyński (2017); Dasenbrook et al. (2015) the non-renewal statistics has a long history in chemical physics where it was used to describe single-molecule processes in spectroscopyCao (2006); Osad’ko and Fedyanin (2011); Budini (2010); Witkoskie and Cao (2006) and kinetics.Saha et al. (2011); Cao and Silbey (2008)

The paper is organised as follows. Section II overviews the derivation of the master equation for electron transport through a vibrating molecular junction. In Section III, we introduce quantum jump operators and derive the expression for the joint probability density of two successive waiting times, $w_{2}(\tau_{1},\tau_{2})$. Section IV describes the results of numerical and analytical calculations. Section V summarises the main results of the paper.

We use natural units for quantum transport throughout the paper: $\hbar=k_{B}=e=1$.

II Master equation in the polaronic regime

The single-molecular junction is a molecule connected to macroscopic source (S) and drain (D) electrodes. The corresponding Hamiltonian is

[TABLE]

The molecule is modelled by Anderson-Holstein model – a single electronic level interacting with a localised vibration. The molecular Hamiltonian is:

[TABLE]

where $\epsilon_{0}$ is molecular orbital energy, $\omega$ is molecular vibration energy, and $\lambda$ is the strength of the electron-vibration coupling. $a^{\dagger}(a)$ creates (annihilates) an electron on molecular orbital, and $b^{+}(b)$ is bosonic creation (annihilation) operator for the molecular vibration. The electronic spin does not play any role here and will not be included explicitly into the equations. Electrodes have noninteracting electrons:

[TABLE]

where $a^{\dagger}_{k\alpha}$ creates an electron in the single-particle state $k$ of the source(drain) electrode $\alpha=S(D)$ and $a_{k\alpha}$ is the corresponding electron annihilation operator. The molecule-electrode coupling is described by tunneling interaction

[TABLE]

where $t_{\alpha}$ is the tunneling matrix element.

Using Born-Markov approximation and Lang-Firsov transormationLang and Firsov (1963) we obtain the master equation Mitra et al. (2004):

[TABLE]

where $P_{nq}(t)$ is the probability that the molecule is occupied by $n$ electrons and $q$ vibrational quanta at time $t$. The transition rates rates are:Mitra et al. (2004)

[TABLE]

– transition from state occupied by one electron and $q$ vibrations to the electronically unoccupied state with $q^{\prime}$ vibrations by the electron transfer from the molecule to $\alpha=S,D$ electrode and

[TABLE]

– transition when electron is transferred from $\alpha$ electrode into the originally empty molecules simultaneously changing the vibrational state from $q$ to $q^{\prime}$. The rates depend on the occupation of electrodes given by Fermi-Dirac numbers

[TABLE]

where $T$ is the temperature and $\mu_{\alpha}$ is the chemical potential of the electrode $\alpha$. The rates also depend on the Franck-Condon factor

[TABLE]

and the electronic level broadening

[TABLE]

where $\rho_{\alpha}$ is density of states in the electrode $\alpha$ taken at molecular orbital energy $\epsilon$.

III Quantum jumps operators for electron tunneling and waiting time distributions

We introduce probability vector ordered in such a way that the electronic probabilities enter in pairs for each vibrational states

[TABLE]

where $N$ is the total number of vibrational states included into the calculations. We also define the identity vector of length $2N$:

[TABLE]

The normalisation of the probability is given by the scalar product between $\mathbf{I}$ and $\mathbf{P}$ vectors $(\mathbf{I},\mathbf{P}(t))$

[TABLE]

Using this probability vector we write the master equation (5,6) in the matrix form

[TABLE]

where $\cal L$ is the Liouvillian operator. The quantum jump operator is $2N\times 2N$ matrix which is defined through the actions on the probability vector:Kosov (2017)

[TABLE]

It describes the tunneling of electron from the molecule to the drain electrode.

We assume that the system has reached the nonequilibrium steady state. Therefore it is described by the steady state density matrix, which is the null vector of the full Liouvillian

[TABLE]

Let us begin to monitor time delays between sequential quantum tunnelings in the nonequilibrium steady state. WTD for two waiting times, $w_{2}(\tau_{2},\tau_{1})$, is defined as joint probability distribution that the first electron waits time $\tau_{1}$ and the next electron waits time $\tau_{2}$ for the tunneling to the drain electrode

[TABLE]

The definition becomes physically obvious if one reads it from right to left: The system is in the steady state described by the probability vector $\mathbf{P}$, then it undergoes quantum jump $J$, then idle without the quantum jump for time $\tau_{1}$, then again undergoes quantum jump $J$, idle for time $\tau_{2}$ and then experiences the quantum jump $J$. WTD for single waiting time between two consecutive tunneling events is

[TABLE]

and again this definition is quite self-explanatory. Let us normalise these distributions

[TABLE]

Here we used that $(\mathbf{I},\;{\cal L}\;\mathbf{X})=0$ for arbitrary vector $\mathbf{X}$. The normalised joint WTD for two waiting time is

[TABLE]

The WTD $w_{1}(\tau)$ has the same normalisation (easy to show by computing the integral over $\tau$)Kosov (2017)

[TABLE]

The formal derivation of WTDs (20) and (21) is presented in appendix A.

Let us also check that these definitions (20) and (21) are consistent with each other. Integrating $w_{2}$ over the second time yields

[TABLE]

Performing integration over the first waiting time (and using ${\cal L}\mathbf{P}=0$) gives

[TABLE]

Therefore, our definitions for single and double time WTDs are consistent with each and have clear probabilistic meaning.

To compute higher-order expectation values and analyse the fluctuations, we introduce the cumulant-generating functions for the joint waiting time probability distribution

[TABLE]

Integrating the cumulant-generating function over $\tau_{1}$ and $\tau_{2}$, we get

[TABLE]

where

[TABLE]

We obtain all possible higher order comulants differentiating $K(x_{1},x_{2})$ with respect to $x_{1}$ and $x_{2}$.

IV Results

IV.1 Equilibrium molecular vibrations

Master equation (5,6) describes the non-equilibrium dynamics of molecular vibrations. Let us first consider the limit where the vibration is maintained in thermodynamic equilibrium at some temperature $T$, which is not necessarily the same as the temperature of electrons in the leads. To implement this limit we use the following separable ansatz for the probabilitiesMitra et al. (2004)

[TABLE]

which assumes that the vibration maintains the equilibrium distribution at all time. The master equation (5,6) is reduced to

[TABLE]

where the vibration averaged rates are defined as

[TABLE]

Let us identify the quantum jump operator for the electron tunneling from the molecule to the drain electrode. We write this jump operator in matrix form and as a dyadic product of two vectors

[TABLE]

Then straightforward vector algebra brings the WTD $w_{1}(\tau)$ (21) to the following form

[TABLE]

For WTD $w_{2}(\tau_{1},\tau_{2})$ (20) we have

[TABLE]

We see that $w_{2}$ is always exactly factorised as a product of two independent $w_{1}$. Therefore, if the molecular vibration is held in thermal equilibrium, then the electronic distribution is always fully reset after each tunneling events and there is no correlation between subsequent tunneling electrons.

IV.2 Nonequilibrium molecular vibrations

Let us now turn our attention to the case of fully nonequilibrium dynamics of molecular vibrations. We must rely on the numerical calculations to get answers in this situation.

We first compute electric current as a function of the applied voltage bias $V_{sd}$. The voltage bias is enforced by shifting symmetrically the chemical potentials of the electrodes $\mu_{S}=V_{sd}/2$ and $\mu_{D}=-V_{sd}/2$. Fig.1 shows the current-voltage characteristics. It has been studied in various details in many works before Galperin et al. (2007); Härtle and Thoss (2011); Cuevas and Scheer (2010) and we show it here simply to serve as a reference - the characteristics steps in the current-voltage characteristic will be related to the behaviour of the waiting time. The steps in the current are due to the resonant excitations of the molecular vibrations by inelastic tunneling of electrons. The steps are observed when the voltage passes through an integer multiple of the vibration energy. We also observe the current suppression in the strong electron-vibration coupling regime due to Franck-Condon blockade.Koch et al. (2006)

Correlations between two subsequent tunneling events will be measured using Pearson correlation coefficient

[TABLE]

Integrals in the Pearson coefficient are computed using cumulant generating function (III): For the correlations we have

[TABLE]

and moments of single waiting time are

[TABLE]

The Pearson correlation coefficient $p$ is widely used in statistics as a measure of correlations between two stochastic variables. It varies between $-1$ and $+1$: $p=0$ means that there is no correlations, $p>0$ indicates positive correlations, and $p<0$ suggests the variables are anticorrelated (negative correlations). In our case positive correlations mean that if the first waiting time increases/decreases then the waiting time for the second electron also increases/decreases. The negative correlations indicate that if the first electron waits longer then the waiting time for the second electron decreases (and vice versa).

Fig.2 shows the Pearson correlation coefficient between successive electron tunneling events (30) computed as a function of the applied voltage bias for different values of electron-vibration coupling. If there is electron-vibration coupling, then the waiting time for successive electron tunneling are not correlated. The correlation does not appear in the weak and moderate electron-vibration coupling regimes (for $\lambda<1$). Only when the electron-vibrational interaction becomes strong $\lambda>1$, the statistical correlations between electronic tunneling events start to emerge.

We focus on the strong coupling regime ($\lambda=2$, $\lambda=3$ and $\lambda=4$ in Fig.2), this is where the voltage dependence of the Pearson correlation coefficient is very interesting. The tunneling events are positively correlated in the narrow voltage window $\omega\leq V_{sd}\leq 2\omega$ (the negative correlations for $\lambda=2$ all have $p<0.02$ – they are statistically negligible). Comparing Fig.1 and Fig.2, we see positive correlations belong to the Franck-Condon blockade regime where the electric current is small and, therefore, the tunneling events are extremely rare.

Let us understand why the electron tunneling events become suddenly positively correlated. First, we compute WTD $w_{1}(\tau_{1})$ for $\lambda=4$ and $V_{sd}=1.56\;\omega$. The dependence of $w_{1}(\tau_{1})$ on the waiting time is shown in Fig.3. The distribution has a large double peaked spike for the fast electrons (the first narrow peak is barely seen on the main figure since it is very close to the axis) - the structure of the spike is zoomed in the insert of Fig.3. To understand the origin of this behaviour of WTD, we compute waiting time distribution for various cutoffs $N$ for the number of the vibrational quanta included in the calculations (Fig.4). For $N=1$ (it means that only states with vibrational quantum numbers $q=0,1$ are included in the calculations), the short time peaks disappear completely. For $N=2$ (vibrational quantum numbers $q=0,1,2$ are included), the short time WTD spike consists only of the one broad peak. For $N=3$ (vibrational quantum numbers $q=0,1,2,3$ are included), the short time WTD spike forms its final two peak structure. Increasing $N$ further does not affect the short time behaviour of WTD at the considered voltage.

We are now in the position to explain the physical reasons behind the appearance of the positive temporal correlations. The mode of the distribution is the first early time peak - it gives us the waiting time for the electrons most often observed in the transport. These fast travelling electrons are responsible for the emergence of the correlations between successive electron tunneling events. We focus on the regime where we observe positive correlations ($\lambda=4$ and $V_{sd}=1.56\;\omega$). After the first quantum jump the first $q=1$ and second $q=2$ vibrational states become populated in the density matrix (57). Then, if the first waiting time $\tau_{1}$ is around the mode of the waiting time distribution $w_{1}(\tau_{1})$ (Fig.3), the vibrational state $q=3$ is predominantly populated after the waiting time $\tau_{1}$(see the appendix for details). If we look at the absolute values of the first 4 components of the Franck-Condon factor $|X_{qq^{\prime}}|,\;q,q^{\prime}=0,..,3$ computed for $\lambda=4$

[TABLE]

we see that it creates the ”shortcut” opening the elastic transport channel through the vibrational state $q=3$. This elastic channel opens only for a short time during the first narrow peak in the WTD $w_{1}(\tau)$ (Fig.3). The higher states $q=4,5,...$ have even larger diagonal Franck-Condon factor but they are not populated on the short time scale at this voltage range. At the higher voltages many vibrational states are populated already in the steady state density matrix, but it does not lead to the correlations since too many elastic channels are available anyway irrespective to the previous tunneling electron.

These temporal correlations between tunneling times should not be confused with ”avalanche” electron transport phenomena in molecular junctions with strong electron-vibrational coupling.Koch and von Oppen (2005); Koch et al. (2006); Lau et al. (2016) The ”avalanche” transport of electrons is observed at intermediate voltage range $3\omega<V_{sd}<6\omega$ just after the Franck-Condon plato in I-V characteristics and it does not involve the correlations between different waiting times.Koch and von Oppen (2005) The avalanche electron transport is simply manifestation of the fact that the mode of WTD is much smaller than the average time given by the same distribution.

V Conclusions

In this paper, we developed the theory for computing joint waiting time distribution for electron transport through a molecular junction with strong electron-vibrational interaction. The molecule is modelled by one molecular orbital coupled with a single localised vibration. We treat electron-vibration interaction exactly and molecule-electrode coupling within the Born-Markov approximation. Using this master equation we computed joint waiting time probability distribution and studied it across various transport regimes to understand the emergence and disappearance of the correlations between successive electron tunneling events.

Our main observations are summarised below:

• •

There are no temporal correlations between subsequent electron tunneling events (the distribution function is completely renewed after each electron tunneling)

1. (a)

for small voltage bias ($<\omega$) and for voltages greater than $2\omega$; 2. (b)

irrespective of voltage bias if electron-vibration coupling $\lambda<1$; 3. (c)

irrespective to voltage and electron-vibration coupling strength if the vibration is maintained in thermodynamic equilibrium.

• •

The temporal correlations between subsequent electron tunnelings emerge, if $\omega<V_{sd}<2\omega$ and $\lambda>1$. The tunneling events become positively correlated which means that the second electron can sense the waiting time for the first electron and if the first electron was fast, then the second electron also would like to be transferred quickly (and opposite, slow to slow correlations are also possible). The observed results are robust, they do not require any special tuning of the parameters other than the physically reasonable choices of applied voltage and electron-vibration coupling.

• •

The physical origin of positively correlated waiting times is the following. After the initial electron transfer to the drain electrode the ground and the first excited vibrational states are predominantly populated. Next, if the first waiting time happens to be around the modal time, then the electron has an opportunity to excite the molecule in $q=3$ state via tunneling through inelastic $q=0$ to $q=3$ or $q=1$ to $q=3$ channels. The excitation of $q=3$ vibrational state creates a ”shortcut” between the source and drain electrodes via the elastic channel with large Franck-Condon factor which can be utilised by the next electron. In other words, if the the first waiting time is the modal time (very short), then the tunneling electron opens $q=3$ elastic channel for swift transfer for the next electron (positive correlations).

Appendix A Waiting time distributions for electron detection events

In this appendix, we define WTDs (21) and (20) using methods of quantum measurement theory. Our derivations follow the theory originally developed in quantum optics to study single photon counting statistics Srinivas and Davies (2010); Zoller et al. (1987); Budini (2010) and extended to quantum transport by Brandes.Brandes (2008)

The Liouvillian in the master equation (15) is decomposed as

[TABLE]

where $J$ is quantum jump operator (16) and ${\cal L}_{0}={\cal L}-J$ generates the evolution of the molecular junction without transferring electrons to the drain electrode. This differential equation is converted to the integral equation

[TABLE]

which then is resolved by iterations

[TABLE]

We introduce operator $M$ - it describes detection of electron transfer from the molecule to the drain electrode (quantum measurement operator). If $\mathbf{P}(t)$ is the probability vector before the measurement, then after the electron detection it becomesBreuer and Petruccione (2002)

[TABLE]

We choose the initial probability vector as

[TABLE]

where $\mathbf{P}$ is the steady state probability vector defined in (17). This choice of the initial state means that we detect the electron transfer to the drain electrode at time $t=0$ in the steady state regime and then we start to monitor the system:

[TABLE]

Using electron detection operator (37), we rewrite (36) in a form which elucidate the probabilistic meanings of its terms:Srinivas and Davies (2010); Zoller et al. (1987); Budini (2010)

[TABLE]

Let us discuss (40). We begin with first term, $e^{{\cal L}_{0}t}M\mathbf{P}$, it is the contribution to the probability vector from all measurements where no electron transfer to the drain electron to occur up to time $t$ after the initial detection at time $t=0$.

The first integral term in this equation can be read as the following: an electron is detected in the drain electrode at time $t=0$ (due to the presence of $M\mathbf{P}$), then no detection of electron is observed up to time $t_{1}$ (due to presence of the ”idle” evolution operator $e^{{\cal L}_{0}t_{1}}$), then the detection of the second electron occurs at time $t_{1}$, and then the system ”idle” without electron transfer to the drain electrode up to time $t$. Therefore, the waiting prefactor $(\mathbf{I},Je^{{\cal L}_{0}t_{1}}J\mathbf{P})/(\mathbf{I},J\mathbf{P})$ must be understood as the probability of observing this process. This waiting factor is exactly our expression (21) for normalized WTD $w_{1}(t_{1})$. The analysis of the second integral term follows exactly the same lines and the waiting prefactor $(\mathbf{I},Je^{{\cal L}_{0}(t_{1}-t_{2})}Je^{{\cal L}_{0}t_{1}}J\mathbf{P})/(\mathbf{I},J\mathbf{P})$ is interpreted as WTD $w_{2}(t_{1}-t_{2},t_{2})$ (20).

Appendix B Probability vectors after quantum jumps

We show in this appendix the probability vectors (normalised) after the quantum jumps associated with electron tunneling to the drain electrode. All probability vectors are computed at the voltage, which corresponds to the maximum Pearson correlation coefficient. The electron-vibration coupling is $\lambda=4$.

We show only first $q=0,1,2,3$ states from the full probability vector (12). Steady state probability and the probability vector (normalized) after the quantum jump are

[TABLE]

Now we compare $Je^{{\cal L}_{0}\tau}J\mathbf{P}$ for different waiting times $\tau=0.05\times 10^{5}$ (corresponds to the first peak in WTD shown on insert plot in Fig.3), $\tau=10^{5}$ (corresponds to the second peak in the WTD shown in insert plot in Fig.3 and $\tau=5\times 10^{6}$ (tail of the WTD):

[TABLE]

From these probability vectors, we see that electron tunneling during the first peak waiting time populates predominantly $q=3$ vibrational state; electron tunneling during the second peak waiting time populates predominantly $q=2$ vibrational state, and then waiting long time to tunnel brings the probability close to the initial $J\mathbf{P}$ vector.