THE JOURNAL OF CHEMICAL PHYSICS 142, 084106 (2015)

Electron transfer statistics and thermal fluctuations in molecular junctions Himangshu Prabal Goswami and Upendra Harbola Department of Inorganic and Physical Chemistry, Indian Institute of Science, Bangalore 560012, India

(Received 6 November 2014; accepted 4 February 2015; published online 24 February 2015) We derive analytical expressions for probability distribution function (PDF) for electron transport in a simple model of quantum junction in presence of thermal fluctuations. Our approach is based on the large deviation theory combined with the generating function method. For large number of electrons transferred, the PDF is found to decay exponentially in the tails with different rates due to applied bias. This asymmetry in the PDF is related to the fluctuation theorem. Statistics of fluctuations are analyzed in terms of the Fano factor. Thermal fluctuations play a quantitative role in determining the statistics of electron transfer; they tend to suppress the average current while enhancing the fluctuations in particle transfer. This gives rise to both bunching and antibunching phenomena as determined by the Fano factor. The thermal fluctuations and shot noise compete with each other and determine the net (effective) statistics of particle transfer. Exact analytical expression is obtained for delay time distribution. The optimal values of the delay time between successive electron transfers can be lowered below the corresponding shot noise values by tuning the thermal effects. C 2015 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4908230]

I. INTRODUCTION

Non-equilibrium quantum transport phenomena is a well known area of research in modern day physics and chemistry. Today, nanofabrication techniques such as molecular beam epitaxy,1,2 lithography,3 and scanning tunneling microscopy4 can be used to control particle transport up to a single atomic level.5 Theoretical models based on several methods such as scattering matrix,6 non-equilibrium Greens function,7,8 quantum master equations,9–11 renormalization groups,12 real-time path integral formalism,13 and random matrix theory14 have been used to analyze the quantum transport. Most of the work in this field has been focused on the average non-equilibrium steady state transport properties of quantum junctions. However, it is known that because of the miniaturization of the system size, quantum interference effects15–17 and thermal fluctuations may influence transport properties18,19 at quantum junctions. The inherent statistical nature of the flux-carrying particles (bosons or fermions) also affects the transport properties at small length scales. Typical experiments measure quantized resistance/conductance using scattering principles and shot noise techniques.20,21 Two probe measurement devices like Hanbury Brown-Twiss or Hong-Ou-Mandel interferometers are used to detect the arrival of particles (photons22,23 and electrons24,25) which demonstrate the effects of inherent quantum nature of particles on the observed statistics. Analysis of photon statistics is based on time average of second order intensity correlation function from which it was inferred that photons tend to arrive in bunches (called bunching) from thermal sources. Electrons are known to exhibit antibunching phenomena,26,27 i.e., in general, they like to avoid one another due to electrostatic interactions. For example, the delay time distribution between two successive electron transfers in a double quantum dot was shown to be totally non Poissonian.28 However, electron bunching has also been 0021-9606/2015/142(8)/084106/8/$30.00

observed in recent experiments with overdense plasma in the attosecond regime,29 free electron lasers,30,31 as well as in electron transport in quantum dot aided by Kondo correlations32 and by shot noise measurements on quantum dots in magnetic field.33 Commonly, bunching and antibunching is a characterization of the statistical analysis of fluctuations in terms of the index of dispersion34 also known as the Fano factor,35–40 based on the first two cumulants of the full probability distribution function (PDF). As a common practice, the analysis of bunched and antibunched statistics in the open quantum systems (sourcesystem-drain approach) are often analyzed under shot noise limit when the electron transfer is unidirectional41,51 and there are no thermal effects. Shot noise corresponds to extreme nonequilibrium limit, i.e., the external bias is infinitely larger than the thermal energy. In a more general case, where the temperature is finite and bias is comparable to the thermal energies, electrons can be transferred between the system and the leads against the applied bias. The unidirectional transfer is affected by the presence of thermal processes and one needs to monitor the net number of electrons transferred to and fro between the system and the leads. The full distribution function of electron transport contains information about the fluctuations to all orders. In general, the analytical form of the full probability distribution function is not known. However, the full distribution function of electron transfer in moletronic devices modeled as single molecule or quantum dot junctions has been obtained numerically to analyze non-equilibrium fluctuations, known as the full counting statistics (FCS).42–50 The FCS of net number of particles transferred is fundamentally connected to the nonequilibrium thermal fluctuations.44,52,53 In the present work, we study effects of thermal fluctuations on the electron transfer statistics and present an analytical result for the corresponding full probability distribution function. For large net number

142, 084106-1

© 2015 AIP Publishing LLC

084106-2

H. P. Goswami and U. Harbola

J. Chem. Phys. 142, 084106 (2015)

(q) of electrons transfer (rare fluctuations), the PDF decays exponentially with q. We analytically identify the rate of decay for the two extrema (q = ±∞) values. The difference in the decay rates at the two extrema is due to the applied bias which breaks the symmetry in the PDF. This asymmetry in PDF is related to the fluctuation theorem. We consider a simple model for electron transfer in an open molecular junction which has been studied earlier.44,53 This consists of two metal leads coupled to a single electronic level. The overall electron transfer statistics in this single resonant level is studied by monitoring the net number of electrons transferred to and fro between the system and the leads. We obtain analytical expressions for the full distribution functions for the net electron transfer statistics by combining the generating function formalism and the large deviation theory. The distribution function is found to fall off exponentially in tails with different rates. The statistics of electron transfer is analyzed in terms of the Fano factor and shown to interchange statistical behavior between bunched, antibunched, and Poissonian by simple manipulation of the thermal fluctuations. We compute the delay time distribution function in presence of thermal effects and show that the optimal delay time can be lowered below the corresponding shot noise value. The paper is organized as follows. In the next section, we introduce the model Hamiltonian and obtain a quantum master equation that describes the dynamics of the system. In Sec. III A, we derive the generating function for the probability distribution of net electron transfer. We present the large deviation theory in Sec. III B, from which we evaluate the probability distribution function for the net electron transfer. The physical parameters responsible for the change in the statistical events of the bunching, random, and antibunching behavior of electrons are discussed in terms of the Fano factor. In Sec. IV, we present analytical form of the delay time distribution and analyze the effects of thermal fluctuations on the optimal delay time. Conclusions are made in Sec. V.

II. FERMIONIC JUNCTION

The quantum junction is made of a single electron level as shown in Fig. 1. This level is connected to two fermionic reservoirs at chemical potentials µl (left lead) and µr (right lead). The Hamiltonian of the system is, Hˆ = Hˆ s + Hˆ l + Hˆ r + Hˆ sν   † ϵ l cˆl†cˆl + ϵ r cˆr†cˆr = ϵ s cˆs cˆs + r

l

+



(Tsν cˆs†cˆν

+ HC),

(1)

sν

where s, r, l, are the system, the right and the left lead orbitals respectively with ν = l,r. Tsν coupling between the system and leads. The operators cˆs†(cˆs ), cˆl†(cˆl ), and cˆr†(cˆr ) represent the electronic creation(annihilation) operators satisfying the Fermi anti commutation relations, { x, ˆ yˆ } = xˆ yˆ + yˆ xˆ such that { cˆi†, cˆj }

= δi j ,

{ cˆi†, cˆ†j }

= { cˆi , cˆj } = 0; i, j = s, l,r.

(2)

The time evolution of the total density matrix, ρˆT (t) is ˆ ρT (t)] and given by the von Neumann equation, ρ˙T (t) = i/~[ H,

FIG. 1. A schematic diagram of a single level system coupled to two Fermionic reservoirs. The chemical potential of the right lead, µ r = µ o . The left lead has higher chemical potential related to the right lead as µ l = µ o + eV . The quantity, eV is the external bias.

A˙ = ∂ A/∂t. Using a second order perturbation in Hˆ sν , invoking the Markovian approximation and the Born approximation (for a generalized derivation see Ref. 10), we obtain a quantum master equation for the reduced system density matrix, ρˆ = Tr B { ρˆT }. This reduction is done assuming the classical nature of the leads which is assumed to have infinite degrees of freedom. The final equation describes the dynamics of the system in the reduced space. We get i ρ(t) ˙ = − [ Hˆ s , ρ(t)] + α cˆs ρ(t)cˆs† − β ρ(t)cˆs cˆs† ~ − α cˆs†cˆs ρ(t) + β cˆs† ρ(t)cˆs + H.C,

(3)

where α = α l + α r and β = βl + βr are the system to leads and leads to system electron transfer rates, respectively, αν = γν (1 − f ν (ϵ s )),

(4)

βν = γν f ν (ϵ s ).

(5)

γν = πnν (ϵ s )|Tsν |2, and nν represent the density of states of lead ¯ ν. f ν (ϵ s ) = (e βν (ϵ s −µν ) + 1)−1 is the Fermi distribution function for the leads and β¯ν = (k BT)−1 is the inverse temperature. For a single electronic level, there are two many-body states, |1⟩ and |0⟩, corresponding to the occupied and unoccupied states of the orbital. There are four pathways of electron transfer (see diagram in Fig. 2). In unidirectional (shot noise) limit, the backward in-rate from the right reservoir, βr , and the backward out-rate from left reservoir, α l , vanish. This corresponds to taking the limit f L (ϵ s ) = 1 and f R (ϵ s ) = 0. In the following, we shall consider all the four processes in our analysis. In the Liouville space, the reduced density matrix is represented by a vector containing the elements | ρ⟩⟩ = { ρ11, ρ00, ρ10, ρ01}, where ρ mm′ = ⟨m| ρ|m ′⟩, m, m ′ = 0, 1. Thus, the quantum master equation in the Liouville space can be represented as, ˆ ρ(t)⟩⟩, where Lˆ is the Liouville operator, | ρ(t)⟩⟩ ˙ = L| −α *. α Lˆ = 2 ... . 0 , 0

β −β 0 0

0 0 +/ 0 0 // . // −iϵ s − (α + β) 0 0 iϵ s − (α + β) -

(6)

084106-3

H. P. Goswami and U. Harbola

J. Chem. Phys. 142, 084106 (2015)

evolve as,





t

| ρ˜(q)(t)⟩⟩ =

tq

dt q 0



×

dt q−1 . . .

0 t2

dt 1Lˆ 2(t q )Lˆ 2(t q−1) . . . Lˆ 2(t 1)| ρ(0)⟩⟩, ˆ

(11)

0

FIG. 2. Ladder diagrams showing the dynamics of the density matrix evolution between the two many-body states |0⟩ (no electron) and |1⟩ (one electron) due to the interaction with the leads. The rates are shown as dotted lines with two arrow heads. There can be four such rates. At time t 1, the operator cˆ †s acts from the left on the bra of the density matrix and cˆ s acts on the ket of the density matrix from the right side. As a result, the density matrix evolves from state |0⟩ to state |1⟩ with the rate, β l . A similar situation is shown at times t 3 and t 5 where the rates are β r . Likewise, the system density matrix jumps from one electron state to no electron state at times t 2 and t 4 with a rate α r when the operator cˆ †s acts from the right on the ket of the density matrix and cˆ s acts on the bra of the density matrix from the left side. The same goes for the time t 6 but at a different rate α l . These processes however happen randomly.

According to Eq. (6), the population dynamics is decoupled from the Fock space coherences, which vanish exponentially with time. We therefore concentrate only on the population dynamics. The new vector, | ρ⟩⟩ consists of the elements, ρ mm only | ρ(t)⟩⟩ ˙ = (Lˆ 1 + Lˆ 2)| ρ(t)⟩⟩,

(7)

where −α Lˆ 1 = 2 * , 0 0 Lˆ 2 = 2 * α ,

0 + , − ββ+ . 0-

(8)

(9)

Here, Lˆ 1 describes the system dynamics in absence of electron transfer and the operator Lˆ 2 induces transitions between different many-body states and creates population changes. In the next section, we use Eq. (7) to investigate the various statistical aspects of non equilibrium electron transport in the quantum junction.

III. ELECTRON TRANSFER STATISTICS A. The generating function

where the time evolution is with respect to Lˆ 1. The physical basis of Eq. (11) can be explained in the diagram shown in Fig. 3. The system initially starts with no electrons until a time t 1 where the off-diagonal Liouvillian, Lˆ 2, acts upon it to transfer an electron between leads and the system. The system then evolves without any transfer of electrons. At t 2, Lˆ 2 again acts to transfer an electron. This process is indefinitely continued till time t q when the q-th electron has been transferred. After that no more electrons are transferred until the observation time t is reached. This is one particular realization that contributes to the process of transferring q electrons. Other processes differ in interactions times t 1, t 2 . . . t q . The net conditional density matrix is obtained by considering all possible values of interaction times t 1, t 2 . . . t q . The probability, P(q, t) of q particle transfer is then obtained by tracing Eq. (11), which can be expressed as  t ˆ Lˆ 2(t q )| ρ˜ q−1(t q )⟩⟩. P(q, t) = dt q ⟨⟨ I| (12) i

0

To evaluate the full counting statistics of electron transfer, we need to close Eq. (12). This is done by defining a moment generating function for the distribution  G(λ, t) = P(q, t)eλq . (13) q

We have identified four elementary processes through which electrons can enter or leave the system (see Fig. 2) and change its state. Let q1(q2) represent the number of electrons that come to (leave) the system from (to) the left lead while q3(q4) represent the number of electrons that come to (leave) the system from (to) the right lead. Then q = q1 + q2 + q3 + q4 in Eq. (13) and the generating function G(λ, t) gives statistics of total number of transfer processes between two leads and the system. However, statistics of the four individual processes can similarly be studied by introducing a vector notation q = {q1, q2, q3, q4} and the four conjugate parameters, λ = {λ1,  λ2, λ3, λ4}, such that λ.q = 4i=1 λi qi . The vector λ represents the counting parameters associated with the transfer of electrons from(to) the system to (from) the leads. Thus, in this

We are interested in the statistics of the net number of electrons transferred between system and reservoirs during time t. We adopt the generating function method as discussed in Refs. 44 and 53. Here, we briefly present this method for completeness. Let | ρ(q)(t)⟩⟩ be the system density matrix such that q electrons are transferred between system and leads, the probability of transferring q particles during a certain time t can be evaluated by tracing this conditional density matrix, ˆ ρˆ(q)(t)⟩⟩, P(q, t) = ⟨⟨ I|

(10)

ˆ is a unit vector. In the interaction picture, the condiwhere | I⟩⟩ tional density matrix, | ρˆ(q)(t)⟩⟩ in Eq. (10) can be shown to

FIG. 3. Feynman diagram of the evolution of the conditional density matrix evolution in the Liouville space under the interaction representation.

084106-4

H. P. Goswami and U. Harbola

J. Chem. Phys. 142, 084106 (2015)

notation, the generating function can be written as  G(λ, t) = P(q, t)eλ.q,

(14)

q

where P(q, t) is obtained from (12) by replacing vector | ρ˜ q−1⟩⟩ with | ρ˜q−I⟩⟩, where I is the unit vector. Combining Eqs. (12) and (14), we get a closed equation for the generating function  ˙ t) = ⟨⟨ I|( ˆ Lˆ 1 + eλi Lˆ 2(i))|G(λ, t)⟩⟩ G(λ, i

ˆ M(λ)|G(λ,t)⟩⟩, ˆ = ⟨⟨ I|

(15)

where Lˆ 2(i) represents the transfer Liouvillian containing the rate  corresponding to the i-th process such that Lˆ 2 = i Lˆ 2(i) and, −α βl eλ1 + βr eλ4 + ˆ M(λ) = 2 * λ2 . (16) λ3 −β , αl e + αr e The eigen values (ζ m ) of this matrix can be written as (m = 1, 2) 
2 ζ m = −(α + β) ± (α + β + 4(αλ βλ − α β), (17) where αλ = α l e + α r e , βλ = βl e + βr e . As λ → 0, ζ1 → 0 and ζ2 → −2(α + β). The formal solution of the generating function in Eq. (15) is given by, λ2

λ3

λ1

electronic transport, this is, however, a difficult task and an exact solution is not possible in terms of simple functions. This problem is, however, not so difficult if we focus only on the steady-state statistics which corresponds to the observation time t being much larger than the relaxation time of the system. For this, we define the following scaled cumulant generating function: 1 (24) S(λ) = lim ln G(λ, t). t→ ∞ t We scale it by 1/t to study the steady-state (time independent) statistics. Substituting the value of G(λ, t) from Eq. (19) in Eq. (24), we get, 1 S(λ) = lim ln[eζ2t ⟨⟨I|R(ζ2)⟩⟩⟨⟨L(ζ2)| ρ(0)⟩⟩ t→ ∞ t + eζ1t ⟨⟨I|R(ζ1)⟩⟩⟨⟨L(ζ1)| ρ(0)⟩⟩]. (25) In the long time limit, the exponential term with the larger eigenvalue (ζ1) dominates the sum inside the logarithm. Thus, in the long time limit, the scaled cumulant generating function is equivalent to the larger eigenvalue, S(λ) = ζ1 = −(α + β) +

λ4

ˆ ˆ M(λ)t | ρ(0)⟩⟩. ˆ G(λ, t) = ⟨⟨ I|e

(18)

We focus only on the net electron transfer between the left lead and the system. We therefore put λ3 = λ4 = 0, λ1 = −λ2 = λ. This allows us to measure the net number of particles q = q1 − q2 transferred between the left lead and the system, irrespective of q3 and q4 values. In the following, we shall therefore drop the vector notations for λ and q. We can re-write the generating function for our system in terms of the left and the right eigen-vectors ⟨⟨L(ζ m )| and |R(ζ m )⟩⟩, of M(λ)  ˆ G(λ, t) = eζ m (λ)t ⟨⟨ I|R(ζ ˆ (19) m )⟩⟩⟨⟨L(ζ m )| ρ(0)⟩⟩,

|R(ζ1)⟩⟩ = ( βl e + βr , α + ζ1), 1 ⟨⟨L(ζ1)| = (α + ζ1, βl eλ + βr ), βλ(ζ2 − ζ1) λ

|R(ζ2)⟩⟩ = ( βl e−λ + βr , α + ζ2), 1 ⟨⟨L(ζ2)| = (−α − ζ1, βl eλ + βr ). βλ(ζ2 − ζ1)

(20) (21) (22) (23)

Note that the generating function in Eq. (19) corresponds to the statistics of net q electrons transferred between the left lead and the system for which we shall compute the FCS in the following. B. Steady-state analysis of fluctuations

The full distribution function, P(q, t) for transferring q electron between lead and the system, in principle, can be directly obtained by inverting Eq. (14), using (19) and then taking the long time limit. Such a direct inversion was done for classical particle transport.54 For the present model of

(α + β)2 + 4 f (λ),

(26)

where f (λ) = α l βr (e−λ − 1) + α r βl (eλ − 1). In order to obtain analytical expression for P(q, t) (here t is much larger than the relaxation time), we make use of the large deviation theory.55,56 The large deviation theory relies on the exponential decay of probabilities of rare events. The PDF is then expressed as,55,56 P( y) ≈ e−tL(y).

(27)

Here, y = q/t is the net number of electrons counted per unit time and hence is the net flux in the system. L( y) is the socalled large deviation function (LDF) defined as a LegendreFenchel transformation of the scaled cumulant generating function L( y) = extλ( yλ − S(λ)).

m

where the left and the right eigen vectors are obtained as,



(28)

The large deviation method is a mathematical approximation where the full PDF, which is the inverse Laplace transform of G(λ, t), is evaluated using the steepest decent method of integration and saddle point approximation.57 As t increases, the large deviation principle’s fidelity keeps increasing and the full distribution can be recovered in the large t limit. Next task is to compute the extremum of the RHS in Eq. (28). In order to solve for the extremum in Eq. (28), we need to solve the following equation: y−

d S(λ) = 0, dλ

(29)

where d 2(α r βl eλ − α l βr e−λ) S(λ) =  . dλ (α + β)2 + 4α r βl (eλ − 1) + 4α l βr (e−λ − 1) (30) Approximating the denominator as  2 f (λ)
(α + β)2 + 4 f (λ) ≈ (α + β) 1 + , (α + β)2

(31)

084106-5

H. P. Goswami and U. Harbola

J. Chem. Phys. 142, 084106 (2015)

we get, d 2(α + β)(α r βl eλ − α l βr e−λ) S(λ) ≈ . dλ (α + β)2 + 2 f (λ)

(32)

Integrating back Eq. (32), we obtain the following approximate expression for S(λ):   2 f (λ) . (33) S(λ) = (α + β) ln 1 + (α + β)2 We shall use this form of S(λ) in the following. Using Eq. (32) and y = q/t, the extremum is found to be at    e F (τ + q)( q˜ + q2 + q)  1  ∗  , (34) λ = ln   2 − q)  2 (τ − q)( q ˜ + q   where F = ln

f r (ϵ s )(1 − f l (ϵ s )) f l (ϵ s )(1 − f r (ϵ s ))

(35)

FIG. 5. The probability distribution function calculated from large deviation theory. The various distributions are simulated at different characteristic times (τ) increasing from left to right. All energy units are in k BTl and eV = 0.5. Γl = Γr = 0.075, E = 0.5, µ o = 0.5, Tr = 1. The exact results are evaluated numerically by inverting the generating function in Eq. (15) and are shown as points. The solid curves represent the distributions from large deviation function in Eq. (37).

is the thermodynamic force coupled to the electron flux and τ 2 − q2 F e , x˜ 2 2 τ x˜ = 2 − (1 + e F ), 2t α r βl

q˜ = 4

(36)

with τ = t(α + β) a dimensionless time. Substituting Eq. (34) in Eq. (28), the LDF is obtained as, qλ∗ − S(λ∗). (37) t The large deviation function in Eq. (37) is shown in Fig. 4 for different bias values as a function of flux y = q/t. The probability distribution function at long times is then obtained by substituting Eq. (37) in Eq. (27). We compare this analytical result for the PDF with the exact numerical result obtained by inverting Eq. (19). This is shown in Fig. 5 for different times, τ. The comparison is quite satisfactory. As we notice, for q = τ, λ∗ diverges. However, the PDF is well defined and decays exponentially with τ, ( )τ 2α r βl P(q = τ,t) ∼ , (38) (α + β)2 L(q, t) =

where the quantity inside the bracket is less than unity and therefore becomes smaller as τ increases. At large q-values (q ≫ τ), the probability distribution decays exponentially: P(q ≫ τ) ∼ e−|q|R , with the rate, √  1 e F ( x˜ + x˜ 2 − 4e F )   if q > 0, ln  √   2 − 4e F 2  x ˜ − x ˜ (39) R= √  −F 2 − 4e F )   1 e ( x ˜ + x ˜   if q < 0.  ln √ 2 x˜ − x˜ 2 − 4e F Thus, the distribution function changes qualitative behavior from a Gaussian, around the maximum, to an exponential in the tails. Note that there is an asymmetry in the tails for the positive and negative q values which is due to the non-equilibrium force F . Equation (39) allows us to compute the ratio of the probabilities in the tails for the positive and negative q values. We obtain (for large |q|), P(q) = eq F . (40) P(−q) This is the manifestation of a deeper symmetry relation which is valid at all q values, q L(q) − L(−q) = − F , (41) t where we have used the symmetry: λ∗(q) + λ∗(−q) = F . The symmetry relation (41) together with the PDF, Eq. (27), leads to the following steady-state fluctuation theorem: lim

t→ ∞

P(q, t) = eq F . P(−q, t)

(42)

Equation (39) shows that even at the equilibrium (F = 0), the probability distribution P(q) decays exponentially for large q. Under the limiting conditions, Tl = Tr , µl = µo + eV and assuming µr = µo = 0, the fluctuation theorem reduces to a more familiar form,44,52,53 FIG. 4. The large deviation function plotted versus number of electrons transferred per unit time, y = q/t. All energy units are scaled by k BTl with E = 0.5, µ o = 0.25. Γl = Γr = 0.075, Tr = 1. The symmetric curve with respect to 0 represents equilibrium (Tl = Tr , eV = 0). From the extreme right, eV values are 3, 2.5, 2, and 0.

qeV P(q, t) = e k BT . (43) P(−q, t) We next analyze the electron transfer statistics in terms of Fano factor. The Fano factor, F, has the information about the

lim

t→ ∞

084106-6

H. P. Goswami and U. Harbola

J. Chem. Phys. 142, 084106 (2015)

dispersiveness of a given data set and describes the inherent statistical nature of events. F > 1, represents a over-dispersed distribution, which corresponds to the positive correlation between events or particles (bunching). Similarly, F < 1 indicates an under dispersed distribution with anti-correlated events (antibunching). F = 1 indicates a totally random (Poissonian) process. The Fano factor is defined as F=

C (2) , C (1)

(44)

where C (i) is the ith order cumulant. Note that the Fano factor as defined above is valid only for processes with C (1) , 0, otherwise it diverges. For the electron transfer statistics, we have ⟨q⟩ = 0 at equilibrium. Thus, we cannot rely on this definition to study equilibrium statistics. In the following, we discuss Fano factor only in the non-equilibrium cases and, based on the Fano factor, show how thermal fluctuations compete with the fluctuations inherent to particle statistics. The i-th order cumulant of P(q, t) is then obtained from the i-th order derivative of S(λ), 

 di (45) C (i) = (q − ⟨q⟩)i = i S(λ) λ=0. dλ Substituting for the first two cumulants and a bit of manipulations, the Fano factor can be expressed as F = 1 − 2ΓE ( f l − f r ) +

2 f r (1 − f l ) , ( fl − fr)

(46)

where ΓE = Γl Γr /(Γl + Γr )2. For f l > f r , as we shall consider throughout this work, all three terms of the left side of Eq. (46) are positive. Clearly, F > 1 (bunching) if ΓE < f r (1 − f l )/( f l − f r )2, F ≤ 1 otherwise, as depicted in Fig. 6. Since ΓE is independent of temperature, the Fano factor can be changed just by manipulating the temperature. Therefore, the statistics of the net number (q) of particles transferred between the system and the lead is strongly influenced by thermal fluctuations. In the unidirectional limit (shot noise regime, T → 0), f l = 1 and

fr = 0, βr and α l vanish and the Fano factor reduces to F = 1 − 2ΓE < 1,

(47)

a well known result which says electron transport through a single resonant level is always antibunched in the shot noise regime. Note that it is only the shot noise limit which describes the inherent statistics of particle transfer. In the case of the net number, q, thermal fluctuations compete with the inherent statistics and may overcome the shot noise which results in the qualitative change in the statistics. The average net current and its fluctuations can be represented in terms of their shot-noise values as (1) C (1) = CSN ( f l − f r ), (2) C (2) = CSN ( fl + fr − 2 fl fr) + 4ΓE2 (Γl + Γr )( f l + f r − f l2 − f r2),

(48)

(1) (2) where CSN = Γl Γr /(Γl + Γr ) and CSN = 2ΓE (Γl2 + Γr2)/(Γl + Γr ) are, respectively, the net current and fluctuations in the shotnoise limit. It is clear that, since 0 ≤ f l , f r ≤ 1, the thermal effects tend to reduce the net current while enhance the fluctuations. This leads to a qualitative change in the statistics. This is, however, not the case with individual particle transfer statistics as we show below. For the four electron transfer processes identified earlier, the Fano factor can be obtained by computing the first two cumulants and using Eq. (44). We notice that the averages and variances of the four processes in general are different. However, there are only three independent process due to the condition that at steady state both the left and the right side fluxes must be the same, i.e., ⟨q1⟩ − ⟨q2⟩ = ⟨q4⟩ − ⟨q3⟩. This relation remains valid for all higher order cumulants. The four Fano factors are obtained as follows:

βα L α βL , F2 = 1 − 2 , 2 (α + β) (α + β)2 α βR βα R F3 = 1 − 2 , F4 = 1 − 2 . (α + β)2 (α + β)2

F1 = 1 − 2

(49)

All four Fano factors are always less than unity, i.e., the individual processes are always antibunched. In the shot-noise limit, the processes q2 and q3 do not exist and F1 and F4 are the same as F given in (47). Thus in the shot noise limit, statistics of all processes is identical while, in general, they are different due to thermal effects. Next, we evaluate analytical expressions for the delay time distribution for electron transport in presence of thermal effects.

IV. DELAY TIME DISTRIBUTION

FIG. 6. The Fano factor, F (Eq. (46)) varied as a function of bias, eV . All units are scaled by k BTl . The various solid curves represent simulations done under varying coupling ratios. From the topmost curve, Γl increases as 0.06, 0.2, 0.5 with fixed Γr = 0.06. µ o = 0, E s = 0.5. The dashed lines represent f (1− f ) the ratio of ( f r− f )2lΓ to ΓE ( f l − f r ) evaluated under the same conditions. E l r Whenever this ratio equals unity, we see F = 1. The inset represents the behavior of the Fano parameter as a function of the coupling ratio. From the right, the curves are evaluated at f l = 0.67, 0.57, 0.47 and fixed f r = 0.3.

Owing to Fermi Dirac statistics, one energy state can be occupied by only one electron. The time spent by an electron in one state and the time difference (delay time) between two successive transfers of electrons to the same state are random quantities and have been studied earlier in the shot noise regime,58 and in mesoscopic conductors in presence of scattering.40 Experimental evidence of such delay time measurements revealed antibunching nature of electrons.28

084106-7

H. P. Goswami and U. Harbola

J. Chem. Phys. 142, 084106 (2015)

However, it is interesting to ask what happens in presence of thermal fluctuations (away from shot noise limit). Suppose, within an observation time window t, an electron is transferred to the system at a time t 1 and departs at time t 2 and then the next electron comes at t 3. The joint probability for the successive arrival of two electrons at times t 1 and then at t 3 is given by  ˆ t L1Lˆ (i)(t 3)Lˆ ( j)(t 2)Lˆ (k)(t 1)| ρ(0)⟩⟩. P(t; t 1,t 2,t 3) = ⟨⟨ I|e ˆ 2 2 2 i, j,k=l,r

(50) At times t 1 and t 3, the transfer Liouvillian has only inward rates, whereas at time t 2, the Liouvillian has the outward rates. Equation (50) can be interpreted as follows. The initial density matrix is subjected to evolution by the transfer of one electron to the system either from the left lead or the right lead at a time t 1 by the transfer Liouvillian containing the respective left or right in-rates ( β). It evolves up to time t 2 without any transfer process. At t 2, the electron is transferred out of the system either from the left or the right lead by the action of the transfer Liouvillian containing the out-rates (α). At time t 3, another electron is transferred to the system by the action of Lˆ 2. The system then evolves up to time t without any transfers. Substituting for Lˆ 1, Lˆ 2 from Eq. (9), the joint distribution (50) is obtained as, P(t; t 1,t 2,t 3) = α β 2 ρ00e−αt e(α−β)t1e(α−β)(t1−t2+t3).

(51)

In order to obtain the delay-time distribution, P(τ++ = t 3 − t 1), where τ++ is the delay time between the arrival of two successive electrons to the system, we need to integrate out t 2(t 1 ≤ t 2 ≤ t 3) and t 1(τ++ ≤ t 1 ≤ t). Finally, integrating out the observation time t from τ++ to ∞ and choosing the steady state value for the initial density matrix, ρ00 = α/(α + β) (See Appendix B), we obtain
αβ e−βτ++ − e−ατ++ , (52) P(τ++) = α− β with the optimal (most likely) value of the delay time given by t ∗ = ln(α/ β)/(α − β).

(53)

The average and the variance of the delay time are ⟨τ++⟩ =

α+ β , αβ

(54)

α2 + β2 . (55) α2 β2 In the shot noise limit, the distribution (52) reduces to the known result58,59
Γl Γr (56) P(τ++) = e−Γr τ++ − e−Γl τ++ , Γl − Γr with the optimal delay-time 2 ⟨τ++ ⟩ − ⟨τ++⟩2 =

ln(Γl /Γr ) . (57) Γl − Γr Comparing (57) with (53), we find that the optimal delay-time ∗ can be reduced further below the shot-noise limit (t ∗ < t SN ) by thermal effects if Γl 1 − fr > . (58) Γr fl ∗ t SN =

FIG. 7. The delay time distribution in the units of k B . The dotted curves represent shot noise and solid curves represent thermally affected distributions. Conditions are generated such that, Eq. (58) is satisfied. A reduction in the maximum delay time is observed. The parameters taken are k B = 1, ϵ s = 0.5, µ 0 = 0, Tl = Tr = 0.5, Γl = 0.22, and Γr = 0.75Γl f l /(1 − f r ). The perpendicular drops indicate numerical evaluation of the optimal delay time.

This is shown in Fig. 7. For the case of symmetric coupling, (Γl = Γr ), the optimal delay-time can be expressed in terms of its value in the shot-noise limit as, ) ( 2 ∗ − 1 ln t f + f +/ , l r (59) t ∗ = SN *. 2 1 − fl − fr , ∗ and therefore for symmetric coupling, t ∗ ≥ t SN , i.e., the thermal effects always increase the optimal delay time.

V. CONCLUSION

We have presented an extensive analysis of the statistics of non equilibrium electronic transport in a quantum junction. An analytical expression for the probability distribution function for the net number (q) of electrons transfer is derived using the large deviation function method. It is interesting to note that the distribution function changes behavior from a Gaussian, which is valid close to the maximum, to an exponential in the tails. The two tails fall off exponentially with different rates. This asymmetry is related to the fluctuation relation and arises due to the non-equilibrium force on the system. The steady state statistics of electron transfer was characterized using the Fano factor. Thermal effects are shown to induce a qualitative change in the statistics by suppressing the average current while enhancing the fluctuations. ACKNOWLEDGMENTS

H.P.G. acknowledges the financial support obtained from the University Grants Commission, New Delhi, India. U.H. acknowledges the financial support from Indian Institute of Science, Bangalore.

APPENDIX A: APPROXIMATION DONE DIRECTLY ON THE CUMULANT GENERATING FUNCTION

In principle, expansion used in Eq. (30) can be done directly on the generating function in Eq. (26). Doing this leads

084106-8

H. P. Goswami and U. Harbola

J. Chem. Phys. 142, 084106 (2015) 11S.

to the following approximate expression for the GF: S(λ) =

2α r βl (e − 1) + 2α l βr (e (α + β) λ

−λ

− 1)

.

(A1)

Using this approximate cumulant generating function, the following expression of the Fano factor is obtained: F = 1 + CE .

(A2)

Thus, F ≥ 1. Clearly, it predicts that antibunching is not possible. This result is not consistent with the exact result in Eq. (46). Making this approximation on the derivative instead, as done in Eq. (33), not only allows us to compute full distribution function but also leads to the exact expression for the Fano factor. APPENDIX B: STEADY STATE DENSITY MATRIX

Here, we evaluate the steady state density matrix of the system which is used in Eq. (51). This is done by solving Eq. (7) with (6). Since the right and left eigen vectors of Lˆ form a complete set, we can evaluate the density matrix by writing  | ρ(t)⟩⟩ = ⟨⟨L(mζ )| ρ(0)⟩⟩|R(mζ )⟩⟩ (B1) mζ

= ⟨⟨L(ζ1)| ρ(0)⟩⟩|R(ζ1)⟩⟩ + e−2(α+β)t ⟨⟨L(ζ2)| ρ(0)⟩⟩|R(ζ2)⟩⟩.

(B2)

Substituting the eigen vectors from Eq. (23) and writing the reduced Liouville space density matrix as | ρ(t)⟩⟩ = { ρ11(t), ρ00(t)}, we get 1 [( β + αe−2t(α+β))ρ00(0) (α + β) + (1 − e−2t(α+β)) β ρ11(0)], 1 ρ00(t) = [(1 − e−2t(α+β))α ρ11(0) α+ β ρ11(t) =

+ (α + βe−2t(α+β))ρ11(0)].

(B3)

Taking t → ∞ limit, we get the steady state values ρ11 = 1A.

β α+ β

ρ00 =

α . α+ β

(B4)

P. Alivisatos, Science 271, 933 (1966). Santori, D. Fattal, J. Vucovic, G. S. Solomon, and Y. Yamamoto, Nature 419, 594 (2002). 3Y. Chen, D. A. A. Ohlberg, L. Xuema, D. R. Stewart, R. S. Williams, J. O. Jeppesen, K. A. Nielsen, J. F. Stoddart, D. L. Olynick, L. Deirdre, and E. Anderson, Appl. Phys. Lett. 82(10), 1610 (2003). 4F. Chen and N. J. Tao, Acc. Chem. Res. 42, 429 (2009). 5M. Fuechsle, J. A. Miwa, S. Mahapatra, H. Ryu, S. Lee, O. Warschkow, L. C. L. Hollenberg, G. Klimeck, and M. Y. Simmons, Nat. Nanotechnol. 7, 242 (2012). 6H. Haug and A. P. Jauho, Quantum Kinetics in Transport and Optics of Semiconductors (Springer, Berlin, 1995). 7S. Datta, Electronic Transport in Mesoscopic Systems (Cambridge University Press, Cambridge, 1997). 8U. Harbola and S. Mukamel, Phys. Rep. 465, 191 (2008). 9U. Harbola and S. Mukamel, Phys. Rev. B 74, 235309 (2006). 10A. Asadian, D. Manzano, M. Tiersch, and H. J. Briegel, Phys. Rev. E 87, 012109 (2013). 2C.

A. Gurvitz and Ya. S. Prager, Phys. Rev. B 53, 15932 (1996). Gezzi, T. Pruschke, and V. Meden, Phys. Rev. B 75, 045324 (2007). 13S. Weiss, J. Eckel, M. Thorwart, and R. Egger, Phys. Rev. B 77, 195316 (2008). 14C. W. J. Beenakker, Rev. Mod. Phys. 69, 731 (1997). 15M. O. Scully and M. S. Zubairy, Quantum Optics (Cambridge University Press, Cambridge, 1997). 16E. Bocquillon, V. Freulon, J. M. Berroir, P. Degiovanni, B. Plaais, A. Cavanna, Y. Jin, and G. Fève, Science 339(6123), 1054 (2013). 17H. P. Goswami and U. Harbola, Phys. Rev. A 88, 013842 (2013). 18R. N. Barnett, C. L. Cleveland, A. Joy, U. Landman, and G. B. Schuster, Science 294, 5542 (2001). 19E. Rossi, J. H. Bardarson, M. S. Fuhrer, and S. Das Sarma, Phys. Rev. Lett. 109, 096801 (2012). 20Ya. M. Blanter and M. Buttiker, Phys. Rep. 336, 1 (2000). 21N. Ubbelohde, C. Fricke, C. Flindt, F. Hohls, and R. J. Haug, Nat. Commun. 3, 612 (2012). 22R. Hanbury Brown and R. Q. Twiss, Nature (London) 177, 27 (1956). 23C. K. Hong, Z. Y. Ou, and L. Mandel, Phys. Rev. Lett. 59, 2044–2046 (1987). 24W. D. Oliver, J. Kim, R. C. Liu, and Y. Yamamoto, Science 284, 5412 (1999). 25H. Kiesel, A. Renz, and F. Hasselbach, Nature 418, 392 (2002). 26M. P. Silverman, Phys. Lett. A 120, 442 (1987). 27M. Iannuzzi, A. Orecchini, F. Sacchetti, P. Facchi, and S. Pascazio, Phys. Rev. Lett. 96, 080402 (2006). 28T. Fuzisawa, T. Hayashi, R. Tomita, and Y. Hirayama, Science 312, 1634 (2006). 29N. Naumova, I. Sokolov, J. Nees, A. Maksimchuk, V. Yanovsky, and G. Mourou, Phys. Rev. Lett. 93, 195003 (2004). 30B. W. J. McNeil and N. R. Thompson, Nat. Photonics 4, 814 (2010). 31G. De Ninno, B. Mahieu, E. Allaria, L. Giannessi, and S. Spampinati, Phys. Rev. Lett. 110, 064801 (2013). 32O. Zarchin, M. Zaffalon, M. Heiblum, D. Mahalu, and V. Umansky, Phys. Rev. B 77, 241303(R) (2008). 33O. Zarchin, Y. C. Chung, M. Heiblum, D. Rohrlich, and V. Umansky, Phys. Rev. Lett. 98, 066801 (2007). 34D. R. Cox and P. A. W. Lewis, The Statistical Analysis of Series of Events (Chapman and Hall, London, 1966). 35U. Fano, Phys. Rev. 72, 26 (1947). 36U. T. Eden and M. A. Kramer, J. Neurosci. Methods 190, 149 (2010). 37G. Schaller, C. Emary, G. Kiesslich, and T. Brandes, Phys. Rev. B 84, 085418 (2011). 38P. Jin, M. Marthaler, J. H. Cole, A. Shnirman, and G. Schon, Phys. Rev. B 84, 035322 (2011). 39M. Kumar, R. Avriller, A. L. Yeyati, and J. M. Van Ruitenbeek, Phys. Rev. Lett. 108, 146602 (2012). 40M. Albert, G. Haack, C. Flindt, and M. Buttiker, Phys. Rev. Lett. 108, 186806 (2012). 41C. Emary, C. Poltl, A. Carmele, J. Kabuss, A. Knorr, and T. Brandes, Phys. Rev. B 85, 165417 (2012). 42D. A. Bagrets and Yu. V. Nazarov, Phys. Rev. B 67, 085316 (2003). 43T. Karzig and F. von Oppen, Phys. Rev. B 81, 045317 (2010). 44M. Esposito, U. Harbola, and S. Mukamel, Phys. Rev. B 75, 155316 (2007). 45A. O. Gogolin and A. Komnik, Phys. Rev. B 73, 195301 (2006). 46G. C. Duarte-Filho and A. M. S. Macedo, Phys. Rev. B 80, 035311 (2009). 47I. Safi and H. Saleur, Phys. Rev. Lett. 93, 126602 (2004). 48D. B. Gutman, Y. Gefen, and A. D. Mirlin, Phys. Rev. Lett. 105, 256802 (2010). 49A. O. Gogolin and A. Komnik, Phys. Rev. Lett. 94, 216601 (2005). 50A. Romito and Y. V. Nazarov, Phys. Rev. B 70, 212509 (2004). 51D. A. Bagrets and Yu. V. Nazarov, Phys. Rev. B 67, 085316 (2006). 52M. Esposito, U. Harbola, and S. Mukamel, Rev. Mod. Phys. 81, 1665 (2009). 53U. Harbola, M. Esposito, and S. Mukamel, Phys. Rev. B 76, 085408 (2007). 54U. Harbola, C. V. den Broeck, and K. Lindenberg, Phys. Rev. E 89, 012141 (2014). 55S. R. S. Varadhan, Comm. Pure Appl. Math. 19, 261 (1966). 56H. Touchette, Phys. Rep. 478, 1 (2009). 57C. M. Bender and S. A. Orszag, Advanced Mathematical Methods for Scientists and Engineers (McGraw-Hill, New York, 1978). 58J. H. Davies, P. Hyldgaard, S. Hershfield, and J. W. Wilkins, Phys. Rev. B 46, 9620 (1992). 59T. Brandes, Ann. Phys. (Berlin) 17(7), 477 (2008). 12R.

Journal of Chemical Physics is copyrighted by AIP Publishing LLC (AIP). Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. For more information, see http://publishing.aip.org/authors/rights-and-permissions.