# Full counting statistics of vibrationally-assisted electronic conduction: transport and fluctuations of the thermoelectric efficiency

## Abstract

We study the statistical properties of charge and energy transport in electron conducting junctions with electron-phonon interactions, specifically, the thermoelectric efficiency and its fluctuations. The system comprises donor and acceptor electronic states, representing a two-site molecule or a double quantum dot system. Electron transfer between metals through the two molecular sites is coupled to a particular vibrational mode which is taken to be either harmonic or anharmonic- a truncated (two-state) spectrum. Considering these models we derive the cumulant generating function in steady state for charge and energy transfer, correct to second-order in the electron-phonon interaction, but exact to all orders in the metal-molecule coupling strength. This is achieved by using the non-equilibrium Green’s function approach (harmonic mode) and a kinetic quantum master equation method (anharmonic mode). From the cumulant generating function we calculate the charge current and its noise and the large deviation function for the thermoelectric efficiency. We demonstrate that at large bias the charge current, differential conductance, and the current noise can identify energetic and structural properties of the junction. We further examine the operation of the junction as a thermoelectric engine and show that while the macroscopic thermoelectric efficiency is indifferent to the nature of the mode (harmonic or anharmonic), efficiency fluctuations do reflect this property.

## I Introduction

Single-molecule junctions offer a versatile playground for probing basic questions in condensed phases physics: How do quantum effects and many-body interactions (electron-electron, electron-phonon) control charge and energy transport processes, thus the operation of nano-scale atomic and molecular devices (1); (2)? How do we accurately and efficiently simulate quantum transport phenomena involving different particles (and quasi-particles), electrons, phonons, spins, polarons? Recent progress in experimental techniques has made it possible to perform sensitive measurements at the molecular scale, in the linear and non-linear transport regimes, to observe signatures of many body effects. Kondo physics, the hallmark of strongly correlated electrons, was observed in different molecules, see e.g. Ref. (3). Coupled electron-vibration processes were probed in single-molecule junctions applying inelastic electron tunneling spectroscopy (4); (5); (6) and Raman spectroscopy tools (7); (8); (9), displaying frequency shifts and mode heating in response to electron conduction. Noise characteristics of the charge current can further expose the nature of the vibrational modes contributing to electron dynamics (10); (11); (5).

Theoretical and computational methodologies dedicated to the effects of electron-phonon interactions on transport in nano-conductors were reviewed in Refs. (12); (13). The celebrated Anderson-Holstein model, with a single electronic orbital coupled to a local phonon mode, exposes an intricate interplay between the electronic and nuclear degrees of freedom. The model has been extensively studied to reveal the behavior of the current and its fluctuations in different regimes of electron-phonon coupling, see for example (14); (15); (16); (17); (18); (19); (20); (21); (22); (78); (24); (25); (26); (27); (28); (29); (30). An extension of the Anderson-Holstein model with a secondary phonon bath was examined in many studies, see e.g. the exploration of thermoelectric transport in a three-terminal junction (31) and the analysis of transient effects (32).

Complementing the Anderson-Holstein model, the donor-acceptor (DA) prototype junction of Fig. 1 allows the exploration of a broad range of problems. The model comprises two electronic states, referred to as “donor” (D) and ”acceptor” (A), following the chemistry literature on electron transfer reactions. Electron transfer between the D and A sites is assisted by a particular vibrational mode, isolated, or coupled to a secondary phonon bath. This model, suggested to describe a molecular electronic rectifier (33), was recently revisited and analyzed using a variety of tools, for example, the Fermi golden rule approach (34); (35), quantum master equations (QME) (36); (37), the non-equilibrium Green’s function (NEGF) technique (38); (39); (40) and from influence functional path integral simulations (41). In molecules, the DA model represents charge transfer between weakly connected chemical groups, facilitated by a vibrational mode (34). In the context of nanoelectromechanical systems (42) the two electronic states can be realized by gate-controlled quantum dots which are coupled by a mechanical system (43), a suspended tunnel-junction such as a carbon nanotube. Electrons transferred between the quantum dots can e.g. excite transverse acoustic modes of the suspended tube (44); (36). The single bosonic mode can also represent a cavity mode assisting electron tunneling between quantum dots, and other hybrid models (45); (46).

The DA model in Fig. 1 offers a rich setting for investigating the role of inelastic-vibrationally assisted electron scattering in far-from-equilibrium (nonlinear) situations. In contrast to coherent conduction, inelastic electron transport can realize nontrivial effects beyond linear response, such as charge and thermal rectification and cross-rectification effects. This was e.g. demonstrated in Ref. (47) within a three-terminal DA junction, by replacing the single vibration by a phonon bath.

Traditionally, quantities of interest in transport experiments and calculations were averaged values for population and currents. However, it should be recognized that nanoscale junctions suffer from strong random fluctuations due to their surrounding environments. It is therefore desirable to develop a probabilistic theory for measurable quantities. Indeed, in small systems the second law of thermodynamics should be replaced by a universal symmetry, the fluctuation theorem (48); (49); (50),

(1) |

where ( is the probability distribution for observing a positive (negative) entropy production during the time interval .

Fluctuations in heat provided to a nanoscale engine and work performed naturally translate to stochastic efficiencies. Universal characteristics of the corresponding probability distribution function , for time-reversal symmetric engines, were explored in Refs. (52); (53); (51), from the principles of classical stochastic thermodynamics. It can be proved that the large deviation function (LDF) for efficiency, defined as , attains a global minimum which coincides with the macroscopic (average) efficiency, and a global maximum which corresponds to the least probable efficiency, coinciding with the Carnot efficiency (52); (53); (51). These universal features are a direct consequence of the fluctuation theorem. Extensions of this analysis to explore efficiency statistics for systems with broken time-reversal symmetry were given in Ref. (54). Fluctuations of the finite-time efficiency were experimentally demonstrated in Ref. (55). Beyond classical thermodynamics, in a recent study the concept of the “stochastic efficiency” was examined within a quantum coherent model of a thermoelectric junction, by employing the non-equilibrium Green’s function technique (56).

In this paper, we provide a complete analysis of charge and energy transport behavior in the donor-acceptor junction by taking a full counting statistics (FCS) approach. We consider two variants of the model as depicted in Fig. 1: Electron transfer may couple to a harmonic vibrational mode (DA-HO model), or to an anharmonic impurity (DA-AH model). The latter case is represented by a two-state system, a truncated vibrational manifold.

We rigorously derive the cumulant generating functions for the models of Fig. 1 under the assumption of weak electron-vibration coupling, handing over the complete information over the models’ steady state transport behavior. From the cumulant generating functions we explore transport in the junctions far from equilibrium, specifically, we aim in identifying transport quantities which are sensitive to the nature of the vibrational mode. Finally, we investigate thermoelectric efficiency fluctuations in the DA junction beyond linear response. Interestingly, our analysis in this paper exemplifies that one can reconcile two central yet disparate techniques, QME and NEGF, to obtain consistent results, within the same order in perturbation theory.

The paper in organized as follows. We introduce the DA junction in Section II and perform a FCS analysis in Sec. III. The cumulant generating function (CGF) of the DA-AH model is derived in Sec. III.1 by employing the quantum master equation approach. The derivation of the CGF for the DA-HO case is detailed in Sec. III.2, using the non-equilibrium Green’s function technique. Two applications are described in Sec. IV: We simulate the junction’s current-voltage characteristics in Sec. IV.1, and examine the statistics of the thermoelectric efficiency in Sec. IV.2, further comparing numerical results far from equilibrium with the linear response (Gaussian) limit. Our findings are summarized in Sec. V. For simplicity, we set throughout derivations.

## Ii Model

The DA junction includes a two-site structure, representing a donor-acceptor molecule (or equivalently, a double-quantum-dot system), placed in between two metal leads. The total Hamiltonian is given by

(2) |

The molecular Hamiltonian includes the donor and acceptor sites,

(3) |

with as a fermionic annihilation (creation) operator at the donor or acceptor sites with energies . The second and third terms in Eq. (2) represent the left () and right () metal leads, modeled by collections of non-interacting electrons,

(4) |

with the fermionic annihilation (creation) operators (. The tunneling energies of electrons from the donor (acceptor) site to the left (right) lead () are included in , and are assumed to be real valued,

(5) |

and represent the Hamiltonians of the molecular vibrational mode and its coupling with the D and A sites (strength ). Assuming an harmonic local mode we write

(6) |

with a bosonic annihilation (creation) operator for a vibrational mode of frequency . The interaction is sometimes refereed to as an “off-diagonal” since electron exchange between the two sites is allowed only via the excitation and/or relaxation of the mode, see Fig. 1. Note that we do not include here a direct tunneling (elastic) term between sites D and A. This simplification allows us to reach closed analytic results for the CGF. The contributing of elastic tunneling processes can be included in an additive manner (47), a reasonable approximation at weak coupling (41).

We diagonalize the electronic Hamiltonian, , and write it down in terms of a new set of fermionic operators ,

(7) |

These operators are related to the original set by (37)

(8) |

with the dimensionless coefficients

(9) |

Here is a positive infinitesimal number introduced to ensure causality. Analogous expressions hold for the set. The expectation values of number operators obey, e.g., at the end, with the Fermi distribution for the lead with chemical potential and an inverse temperature . Using the new operators, we write down the total Hamiltonian for the DA-HO junction as

(10) | |||||

The last term in the Hamiltonian describes electron-hole pair generation assisted by the interaction with the vibrational mode.

In a simpler version of this model we replace the infinite level spectrum of the harmonic oscillator by a truncated two-level system, to mimic a highly anharmonic vibrational mode. The Hamiltonian of this DA-AH model can be conveniently written in terms of the Pauli matrices as

(11) |

The DA-HO and DA-AH models described above are simple enough to allow us to derive expressions for the corresponding CGFs. Meanwhile, we (i) reconcile the NEGF approach with QME, to assist in method development in the area of quantum transport, and (ii) provide analytic results for transport characteristics, to bring intuitive guidelines for functionality.

We substantiate our modeling by describing connections to ab-initio studies of molecular conduction (57). The donor and acceptor sites could represent spatially separated electronic states in a molecule or atoms in an atomic wire, with the electron-phonon interaction of Eq. (6) describing a bond-length stretching mode. The rigid motion of a molecule/atomic chain between the metal leads can be modeled by an additional electron-phonon interaction Hamiltonian . This mode is not included in the present study. As well, we ignore direct tunneling between electronic sites in the form . Detailed ab-initio studies of electron-phonon inelastic effects in molecular junctions prepare direct tunneling elements, frequencies of active modes, and their electron-phonon matrix elements, see e.g. Refs. (58); (57); (59); (60); (61). Certain type of molecular junctions could be well represented by our model, when electron transport due to dominates over both elastic tunneling and the diagonal electron-phonon coupling . This is the case e.g. in Ref. (62), considering charge transfer through a biphenyl molecule with the torsion motion assisting electron hopping between the (almost orthogonal) benzene rings. Particularly, when , calculations of transport in double quantum dot systems show that the inelastic component of the current can dominate the elastic term (47).

We also justify our model in the language of molecular orbitals, electronic eigenstates of the molecule. Consider two orbitals, each coupled to both metal leads- but in an asymmetric manner: one orbital couples strongly to the left lead but weakly to the right, the other molecular orbital is strongly coupled to the right lead but weakly to the left side. In the absence of electron-phonon interaction this molecule supports very small currents. It will however turn into a good conductor at high enough temperatures when phonons contributing to (6) are active, supporting inelastic current. For an extended discussion, see Ref. (63).

Turning to the the AH model, the two-state impurity describe deviations from the harmonic picture. Besides electron-phonon coupled situations, the model could represent electron transport junctions interacting with a local spin impurity, see e.g. Refs. (64); (65), demonstrating electronic read-out of nuclear spins.

We finally comment that in the non-crossing approximation, when elastic and inelastic tunnelling events do not interfere (66); (39), the current (or more generally, the cumulant generating function) can be written as a sum of elastic and inelastic terms. While here we treat the inelastic component only, the CGF for elastic transport is well known (67), bringing in the standard Landauer formula for currents.

In the next Section we develop a counting statistics approach for charge and energy transfer processes. We then analyze first the (simpler) DA-AH model using a QME approach, followed by the investigation of the DA-HO junction by utilizing the NEGF technique.

## Iii Counting statistics for charge and energy currents

The cumulant generating function contains information over the statistics of transferred particles and energy flowing across the system, potentially far from equilibrium. Here we are interested in the CGF for coupled particle (charge) and energy currents. Such a two-parameter CGF is necessary for obtaining later in Sec. IV.2 the statistics of efficiency in a thermoelectric engine.

We define the particle () and energy () current operators from the rate of change of electron number and electron energy in one of the leads, say , and write

(12) |

Here is the number operator for the total charge in the right compartment (right lead plus attached acceptor site). Similarly, . The operators are written in the Heisenberg () picture, thus they should be evolved with the total Hamiltonian for either model, where . We follow the convention that the current flowing out of the right lead is positive. Changes in the total energy and electron number in the lead during the time interval ( and are initial and final observation time, respectively), are given by the integrated currents.

(13) |

Since at any instant , it is possible to construct the so-called “characteristic function” , corresponding to the joint probability distribution for the charge and energy currents. Following the two-time measurement procedure one can define the characteristic function as (49); (50); (68); (69)

(14) |

where and are the counting fields for energy and particles, respectively. represents an average with respect to the total density matrix at the initial time. We assume that , a factorized-product form for the electronic degrees of freedom and for the vibrational part. The leads are maintained in equilibrium at temperature and chemical potential , , and the states are described by the grand canonical distribution function , with as the grand canonical partition function. Eq. (14) can be reorganized as

(15) | |||||

The second line introduces the definition of the total, counting-field dependent, density operator. We trace out its electronic degrees of freedom, () and express the characteristic function in terms of the reduced density matrix for the vibrational mode,

Note that the forward and backward evolution operators are not hermitian conjugates. For example, the forward propagator is

(17) |

with the counting-field dependent total Hamiltonian

is a system operator; the model is reached when , . The model is realized with and . The electron-phonon coupling term here depends on the counting field. Herein, we use as a short-hand notation for both and . To facilitate our discussion below, we define the operators as

(19) |

These operators correspond to the bath operator coupled to the system, see Eqs. (10) and (11), now dressed by counting fields as a consequence of the measurements of charge and energy. Note that the sign convention for corresponds to the respective time evolution operator. The counting-field Hamiltonians [Eq. (LABEL:eq:Hlam) and the complementarity term for the backward evolution] can be organized in a form convenient for a perturbation expansion in ,

(20) |

with incorporating the two metals with the hybridized states, see Eq. (7).

### iii.1 DA-AH model: Quantum master equation approach

We derive the counting field dependent quantum master equation (49); (71) for the model (11) under the assumption that the coupling between electron-hole pair generation and the vibrational mode is weak (72). Unlike other standard QME approaches for molecular junctions, in which the molecule-metal coupling is considered weak (73); (56), in the present derivation the metal-molecule hybridization [defined below Eq. (27)] can be made arbitrarily large, absorbed into the leads spectral density by the exact diagonalization procedure presented in Sec. II. Taking the time-derivative of Eq. (LABEL:eq:vib) we get

(21) |

The operators here are written in the interaction representation, , with including the uncoupled electrons and vibration. By formally integrating this equation we receive the exact form (),

(22) | |||||

We now follow standard steps as in the derivation of the weak coupling-Markov Redfield equation. The initial condition is assumed fully factorized, is replaced by the initial condition, , and the upper limit of integration in extended to infinity, assuming Markovianity of the electron baths. The equation of motion for the counting field dependent reduced density matrix (describing the dynamics of the vibrational mode) depends on the following relaxation () and excitation () rates

(23) | |||||

Here . An explicit calculation of the relaxation rate gives

(24) | |||||

The first term represents an inelastic process with an electron hoping from the left lead to the right one, by absorbing one quanta , satisfying energy conservation with . This process, which goes against our convention of a positive current (flowing right to left), contributes negative charge and energy currents, reflected by the negative sign for and in the exponent. The second term in Eq. (24) corresponds to the reversed process with electron hopping from the right metal to the left, observing , with a positive sign for and . This downward rate assists in cooling the vibrational mode. In the complementary excitation process electrons lose energy to the vibration, heating up the junction. Eq. (24) can be decomposed into two separate contributions,

(25) |

and similarly for . We define the spectral densities for the baths (metals) as

(26) |

Using the transformations (9) we note that these functions are Lorentzian-shaped, centered around the donor ( or acceptor () site energies,

(27) |

with . In terms of these spectral functions, bath induced relaxation rates (24) are given by

(28) |

We also define the rates from Eqs. (23)-(28), only missing the identifier. The rates are nonzero as long as: (i) for a given energy, both left and right leads are not fully occupied or empty and (ii) the overlap between the spectral densities differing by one quanta of energy is non-negligible. Note that because of the weak electron-phonon coupling approximation, each electron tunnelling process involves absorption/mission of a single quanta of energy. In other words, the dynamics is completely described by single-phonon excitation and relaxation rates.

As mentioned above, we apply the weak-coupling Born Markov approximation on Eq. (22). By further ignoring off-diagonal coherence elements for the reduced density matrix, we obtain the population dynamics for the vibrational states (written here for an arbitrary number of levels, ),

(29) | |||||

where and denotes the -th vibrational level. and are rates evaluated at . For the DA-AH model, , it can be shown that off-diagonal coherences do not appear in the Born-Markov approximation without further assumptions (37), and the dynamics of the population follows

(30) |

These equations can be written in a matrix form as

(31) |

where . The long-time (steady state) limit defined as

(32) |

provides the CGF, where is the identity vector. In this limit, only the smallest eigenvalue of the Liouvillian survives. The final result for the CGF for the DA-AH junction is given by

We label this CGF by ‘AH’, to highlight the mode anharmonicity. Recall that collects two counting fields and , for charge and energy, respectively. The CGF satisfies the fluctuation symmetry

which can be verified by examining the rates in Eq. (24). Under the transformations and the rates transform as

(35) |

The extra factors cancel out in Eq. (LABEL:eq:CGF-AH), to satisfy the fluctuation symmetry.

The charge and energy currents and the corresponding zero frequency noise powers can be readily obtained by taking derivatives of the CGF with respect to the counting fields. For example, the particle () and energy () currents are obtained from

(36) |

The zero frequency noise of these currents are

(37) |

where is the second-cumulant.

### iii.2 DA-HO model: Non-equilibrium Green’s function approach

We now focus our attention to the DA-HO junction of Fig. 1(b), described by Eq. (10). In this case, electron-hole pair excitation is coupled to a harmonic oscillator mode. We employ the NEGF technique (74); (75) to derive the CGF of this model, again correct up to the second order of the electron-phonon coupling parameter . Note that the presence of an infinite number of vibrational levels for the HO mode makes it difficult to obtain a closed form for the CGF under the QME approach, since the Liouvillian is an infinite-dimensional matrix.

We begin with Eq. (15), now identifying the initial time by , and write down the characteristic function on the Keldysh contour (see Fig. 2) as (69)

(38) | |||||

The forward upper (backward lower) branch of the Keldysh contour corresponds to the modified unitary evolution (). Evolving the counting fields on two branches of the Keldysh contour with two different signs is the main essence of counting statistics problems. The normalization condition is trivially satisfied with . In the above expression is the contour ordered operator, which orders operators according to their contour time argument; earliest-time operators appear at the right. is the contour time dependent counting function. In the upper branch, , in the lower branch, .

Moving to the interaction picture with respect to the Hamiltonian , we write the characteristic function as

(39) |

where the counting field dependent interaction term is

introducing the short notation . Our objective is to calculate the CGF, correct up to second-order in electron-phonon coupling but exact to all orders in the metal-molecule hybridization. It can be shown that a naive perturbative calculation of Eq. (39) in terms of the electron-phonon coupling leads to a violation of the non-equilibrium fluctuation symmetry. Moreover, such a perturbative treatment cannot capture the correct non-equilibrium phonon distribution, and it will lead to (an incorrect) long-time solution for the vibrational density matrix which depends on the initial- arbitrary state .

In order to restore the fluctuation symmetry, and obtain the correct non-equilibrium phonon distribution [which depends on the temperatures of the electronic baths and their chemical potentials, but not on ], one has to sum over an infinite subclass of diagrams in this perturbative expansion. This procedure takes into account all electron scattering processes which are facilitated by the absorption or emission of a single quanta . Physically, the summation collects not only sequentially tunneling electrons, but all coordinated multi-tunneling processes, albeit with each electron interacting with the mode to the lowest order, to absorb/emit each a single quanta . This summation can be achieved by exploiting the random-phase approximation (RPA) as done in Refs. (27); (76); (77), also referred to as the self-consistent Born approximation (57). Summing over a particular set of diagrams (ring type) in the perturbative series, see Fig. 3, we reach the following expression,

The symbol denotes an integration over contour time variables . For example,

(41) |

Here, is the identity matrix in the Keldysh space and is the free phonon Green’s function,

(42) |

with , proportional to the phonon displacement operator. is the counting-field dependent electron-hole Green’s function. It describes electron hopping processes from the left to the right lead, and vice versa,

This expression is symmetric under the exchange of the contour time parameters and . Recall that the tilde symbol advices that the Green’s function is dependent. This propagator involves free electron Green’s functions for the left and right leads,

(44) |

Explicit expressions for different components of these Green’s functions in real time are given in Appendix A. Here, we write down the lesser component, given as,

(45) | |||||

The greater component is obtained from . It is clear that the greater (lesser) component corresponds to the electronic bath-induced transition rates within the vibrational mode, (), see definitions (23) and a physical explanation below Eq. (24).

Projecting Eq. (LABEL:CGF-RPA) to the real time and invoking the steady state limit by taking , we write the CGF as (75); (68),

(46) | |||||

where with as the third Pauli matrix. Note that we renormalized the CGF with a counting field independent term . The matrix can be written explicitly as

where stand for the retarded, advanced, time-ordered, anti-time ordered, lesser and greater components of the Green’s functions, respectively. The free-phonon Green’s functions