# Vibrational cooling, heating, and instability in molecular conducting junctions: Full counting statistics analysis

## Abstract

We study current-induced vibrational cooling, heating, and instability in a donor-acceptor rectifying molecular junction using a full counting statistics approach. In our model, electron-hole pair excitations are coupled to a given molecular vibrational mode which is either harmonic or highly anharmonic. This mode may be further coupled to a dissipative thermal environment. Adopting a master equation approach, we confirm the charge and heat exchange fluctuation theorem in the steady-state limit, for both harmonic and anharmonic models. Using simple analytical expressions, we calculate the charge current and several measures for the mode effective temperature. At low bias, we observe the effect of bias-induced cooling of the vibrational mode. At higher bias, the mode effective temperature is higher than the environmental temperature, yet the junction is stable. Beyond that, once the vibrational mode (bias-induced) excitation rate overcomes its relaxation rate, instability occurs. We identify regimes of instability as a function of voltage bias and coupling to an additional phononic thermal bath. Interestingly, we observe a reentrant behavior where an unstable junction can properly behave at a high enough bias. The mechanism for this behavior is discussed.

## I Introduction

Can molecules serve as reliable components in electronic circuits? A major obstacle in realizing molecular-based electronic devices is junction heating and breakdown, the result of vibrational excitation by the electron current (1); (2); (3); (4); (5); (6); (7); (8); (9); (10). This situation generally occurs once the bias voltage exceeds typical molecular vibrational frequencies and the electronic levels are situated within the bias window. If energy dissipation from the conducting object to its environment (metals, solvent) is not efficient, the molecular conductor experiences significant heating, ultimately leading to junction breakdown. A related question, the possibility for a nonequilibrium induced cooling of the junction has been the topic of recent experimental and theoretical studies (11); (4); (12); (13).

In this paper, we study the problem of bias-induced molecular cooling, heating, and (potential) junction breakdown due to vibrational instabilities, using the Donor (D)-Acceptor (A) Aviram-Ratner electronic rectifier setup (14), see Fig. 1. By coupling electronic transitions within the junction to a particular internal molecular vibrational mode, significant molecular heating can take place once the donor level is lifted above the acceptor level, as the excess electronic energy is used to excite the vibrational mode. This process may ultimately lead to junction instabilities and breakdown (16). The model can also demonstrate current-induced cooling at low bias, when tuning the junction’s parameters.

Within this simple system, several issues are of fundamental and practical interest. First, one would like to understand the role of mode anharmonicity in the transport process and in the heating or cooling behavior. Second, the molecular vibration under investigation, the one controlling junction stability, can be assumed to be well isolated from other modes. Alternatively, this mode may be coupled to other phonons, allowing for energy damping to a larger environment. These two situations should result in distinctive cooling or heating behaviors. These issues will be explored here. Other relevant challenges which are not considered here are the possibility to selectively excite vibrational modes in the molecule, using voltage bias (17), or more generally, to drive molecular motion or trigger chemical dynamics (18).

Using a full counting statistics (FCS) approach, our analysis further contributes to the institution of fluctuation relations in open many-body quantum systems. Fluctuation theorems (FT) for entropy production quantifies the probability of negative entropy generation, measuring “second law violation” (19); (20). Such “anomalous” processes are relevant at the nanoscale. While originally demonstrated in classical systems (21), recent experimental efforts are dedicated to explore their validity within quantum systems (22). From the theoretical side, the extension of the work and heat FT to the quantum domain has recently attracted significant attention (23); (24). Specifically, a quantum exchange FT, for the transfer of charge and energy between two reservoirs maintained at different chemical potentials and temperatures, has been derived in Ref. (25) using projective measurements, and in Ref. (23) based on the unraveling of the quantum master equation. It is of interest to test these relations in particular cases, e.g., for systems strongly coupled to multiple reservoirs, when the reservoirs cooperatively affect the subsystem (26), including nonmarkovian reservoirs (27); (28); (29), and in models showing coupled charge and energy transfer processes, yet the respective fluxes are not tightly coupled. The system investigated here corresponds to the latter case.

Different flavors of the phonon-assisted-tunneling model have been analyzed in the literature (9). Among the various techniques adopted we list solution of the dynamics as a scattering problem (30), extension of the basic nonequilibrium Green’s function formalism to include molecular vibrations (31), or the use of master equation approaches (32). In this paper, we exploit the latter method, and present a full-counting statistics of the system, allowing for the exploration of charge current, energy current and noise processes at the same footing. Further, we analytically obtain the cumulant generating function (CGF) of the model, allowing for the verification of the steady-state charge-energy fluctuation theorem in this many-body quantum system.

The objectives of this work are therefore twofold: (i) to analyze a simple model that can elucidate cooling, heating and instability mechanisms in molecular rectifiers, specifically, to understand the roles of mode anharmonicity and additional damping routes, and (ii) to establish the steady-state entropy production fluctuation theorem within a nonequilibrium quantum model, transferring charge and energy between the reservoirs in a cooperative manner. Recent studies have analyzed the role of electron-vibration interaction on the full counting statistics (FCS) within different approaches (33); (34); (28); (35); (36). Complementing these efforts, our treatment offers an analytic structure for the CGF, allowing for a clear inspection of the microscopic processes involved.

The plan of the paper is as follows. In Sec. II we introduce the D-A molecular rectifier and its two flavors, either including a harmonic or an anharmonic internal vibration. In Sec. III the anharmonic model is analyzed within a FCS approach, demonstrating cooling, heating and instability dynamics at different parameter regions. The case with an additional phonon bath is considered in Appendix A. Sec. IV explores the harmonic mode model. Sec. V concludes.

## Ii Model

Our model includes a biased molecular electronic junction and a selected internal vibrational mode which is coupled to an electronic transition in the junction. This mode possibly interacts with other (reservoir) phonons, an extension presented in Appendix A. For a schematic representation, see Fig. 1. Generally, this model allows one to investigate the exchange of electronic energy with molecular (vibrational) heating. The total Hamiltonian is given by the following terms,

(1) |

The first term, , stands for the molecular electronic part including two electronic states

(2) |

Here, () is a fermionic creation (annihilation) operator of an electron on the donor or acceptor sites, of energies . The second and third terms in Eq. (1) describe the two metals, , , each including a collection of noninteracting electrons

(3) |

The hybridization of the donor state to the left () bath, and similarly, the coupling of the acceptor site to the right () metal, are incorporated into ,

(4) |

The Hamiltonian further includes an internal molecular vibrational mode of frequency . The mode displacement from equilibrium is coupled to an electron hopping in the system with an energy cost , resulting in heating and/or cooling effects,

(5) |

Here, () represents a bosonic creation (annihilation) operator. Note that in our construction the donor and acceptor sites are coupled to each other only through the interaction with the vibrational mode. We now diagonalize the electronic part of the Hamiltonian, , to obtain, separately, the exact eigenstates for the -half and -half of ,

(6) |

Assuming that the reservoirs are dense, their new operators are assigned energies same as those before diagonalization. The donor and acceptor (new) energies are assumed to be placed within the band of continuous states, excluding the existence of bound states. The old operators are related to the exact eigenstates by (37)

(7) |

where the coefficients, e.g., for the set, are given by

(8) |

Similar expressions hold for the set. It is easy to derive the following relation,

(9) |

with the hybridization strength . The expectation values of the exact eigenstates are

(10) |

where denotes the Fermi distribution function. An analogous expression holds for . The reservoirs temperatures are denoted by ; the chemical potentials are . With the new operators, the Hamiltonian (1) can be rewritten as

(11) | |||||

In this form, the model generally describes the process of an electron-hole pair excitation by a molecular vibration. We denote it by , to highlight the vibrational mode harmonicity. A simple version of the model is reached by replacing the harmonic mode by a two-state system (spin), using the Pauli matrices,

(12) | |||||

The truncated harmonic spectrum imitates an anharmonic mode, as only several (two in the present extreme case) states are bounded within the anharmonic potential (38). We denote this Hamiltonian by , to indicate on the anharmonicity of the molecular mode. The dynamics of this model should coincide with the behavior dictated by , at low temperatures.

Charge and energy transfer dynamics in these models can be followed by studying electronic properties (39); (35); (36). In contrast, here we explore the junction response to an applied voltage bias by studying the vibrational mode excitation and relaxation dynamics. The analysis of the two-state model, Eq. (12), therefore turns out to be simpler than the case when the vibrational mode has an infinite spectrum. In what follows, we derive in details the CGF for the anharmonic-mode case. Appendix A generalizes this calculation to include an additional dissipative thermal bath. We then extend these results and discuss the model conveyed by Eq. (11).

## Iii Anharmonic-mode Rectifier

### iii.1 Impurity dynamics

We explore the dynamics of an anharmonic mode, referred to as an “impurity”, or a two-state-system (TLS), within an electronic rectifier, assuming a weak donor-acceptor - mode interaction. We rewrite Eq. (12) as

(13) |

by defining the electron-hole pair generation operator as

(14) |

The Hamiltonian (13) can be transformed to the spin-fermion model (40); (41) of zero energy spacing, using the unitary transformation

(15) |

with . The transformed Hamiltonian is given by

(16) |

In this form, the TLS dynamics can be simulated exactly using an influence-functional path integral approach (42).

Back to (13), we denote the TLS ground state and excited state by and , with energies and , respectively. We express the Pauli operators by these states, , . Next, using the quantum Liouville equation, we obtain kinetic-rate equations for the states population (=0,1) (43); (38). This standard derivation involves a second order perturbation theory treatment with respect to , the mode-molecule coupling parameter, followed by a Markov approximation. The resulting equation for the reduce density matrix take the simple form ()

(17) | |||||

Here, represents the (mode-molecule) coupling term in Eq. (13). The operators are written in the interaction representation, and we trace over the electronic degrees of freedom. The reservoirs are maintained in a grand canonical state as ; is the partition function of the bath. Identifying the diagonal matrix elements as population, , we obtain the kinetic equation

(18) |

In this model, the off-diagonal elements of the reduced density matrix are naturally decoupled from the population dynamics (44). The excitation () and relaxation () rate constants are given by Fourier transforms of bath correlation functions

(19) |

enclosing electron-hole pair excitation processes,

The operators are given in the interaction representation, e.g., . As we separately trace over the and -baths’ degrees of freedom, it can be shown that the rate constants can be decomposed into two contributions,

(21) |

satisfying

Similar relations hold for the right-to-left going excitations. The energy in the Fermi function is measured with respect to the (equilibrium) Fermi energy, placed at , and we assume that the bias is applied symmetrically, . The four rate constants describe distinct electron-hole excitation processes, depicted in Fig. 2. At forward bias, if we set the effective density of states (DOS) of the bath to lie higher in energy that the DOS of the right bath, we immediately note that the rate should dominate over , potentially leading to ”population inversion” of the vibrational mode. Utilizing electronic reservoirs with energy dependent DOS is thus the basic ingredient of the instability formation here, as we show below. For convenience, we define the spectral density for the bath as

(23) |

Explicitly, using Eq. (8), we find that this function has a Lorentzian lineshape, and that it is centered around either the D or A level,

(24) |

Using the spectral density function, we express the terms in Eq. (LABEL:eq:rate1) as integrals

The following relations hold ( and ),

(26) |

In equilibrium, detailed bias is therefore maintained.

The dynamics conveyed by Eqs. (18)-(LABEL:eq:rate2) is non-separable in terms of the two metals, in contrast to simple linear interaction cases (38). In other words, the reservoirs cooperatively excite or damp energy from the impurity, thus their action is non-additive.

It should be noted that while we assume a weak interaction limit, between electron-hole pair generation and the vibrational mode, our scheme does not enforce weak metal-molecule coupling; this part is exactly diagonalized to yield the reservoirs spectral function, peaked about the D or A levels. If one where to force weak metal-molecule interaction, the spectral functions (24) would reduce to delta functions, and , and the resulting rates would be evaluated at the donor and acceptor levels, e.g., This also implies that charge and energy currents are not “tightly coupled” here, such that for each transferred electron not necessarily precisely one quanta of energy should be gained or drained at either contact. In this aspect, our study complements the work reported in (10). There, using the small polaron transformation, the coupling of the molecular bridge to the leads is assumed to be weak, while its coupling to the vibrational mode can be made large.

### iii.2 Resolved charge and energy equations

We write here a closed expression for the cumulant generating function, following the approach developed in Refs. (26); (27). It will allow us to obtain the current, its noise power, and to confirm the FTs in this system. We define as the probability that by the time the impurity (TLS) occupies the state , electrons have been transferred from the metal to the side, and a net energy has been transferred, to . Resolving Eq. (18) to its charge and energy components, we find that this probability satisfies the following equation of motion (26); (27),

(27) | |||||

One could reason this rate equation as follows. In the first equation, the term stands for the decay rate of . The second line describes a process where by the time the TLS occupies the ground state, excess electrons have arrived at the terminal, and an overall of energy has been absorbed at the bath. At the time an electron-hole pair excitation generates an electron at the bath, leaving a hole at the metal. This charge transfer process is accompanied by an electronic energy transmission at the amount of : An electron leaving the bath has a total energy , however only is gained by the bath. The rest, at the amount of , is gained by the vibrational mode. A similar reasoning can explain other terms in Eq. (27). For convenience, the factor in Eq. (LABEL:eq:rate2) has been absorbed into the definition of .

We Fourier transform the above system of equations with respect to both charge and energy, to obtain the characteristic function . It depends on the energy counting field and the charge counting field ,

It satisfies the differential equation

(29) |

where the matrix contains the following elements

Here,

The cumulant generating function is formally defined as

(32) |

where we introduced the short notation , that is the probability to transfer by the time , electrons and an energy from left to right, irrespective of the state of the TLS. The charge and heat currents can be readily derived, by taking the first derivative of the CGF with respect to either or ,

(33) |

The quantity denotes the total energy transferred from to by the (infinitely long) time ; similarly counts the particles (electrons) transferred in the same direction, by that time. The zero frequency noise current power density can be similarly obtained,

(34) |

The CGF can be expressed in terms of as

(35) |

with , a left vector of unity. It is practically given by the negative of the smallest eigenvalue of the matrix ,

are the matrix elements of , see Eq. (LABEL:eq:mu).

### iii.3 Fluctuation theorem

We confirm next the following symmetry

(37) |

with . In order to prove this, we focus on the product in Eq. (LABEL:eq:QE),

(38) | |||||

Under the transformation and , using the relation , we find that

(39) | |||||

Similarly, one could show that

(40) |

The extra factors cancel, and we recover the symmetry

(41) |

confirming Eq. (37). We can now demonstrate the validity of a fluctuation relation for this non-equilibrium system. The probability to transfer the energy by the long time , from to , is given by the inverse Fourier transform of Eq. (32),

with . Similarly, the quantity represents the probability that charged particles and an energy have been transmitted in the opposite direction, right to left, up to time . Based on the symmetry Eq. (37), one can show that (23)

(43) |

which is often written in a compact form as

(44) |

This expression goes beyond standard metal-molecule weak-coupling schemes as the energy and charge transfer and not tightly coupled, and the energy can take continuous values, unlike Refs. (46); (47); (39).

It should be noted that the above derivation has assumed charge and energy conservation between the two reservoirs. The full particle-energy counting statistics, without such an assumption, would begin with the probability distribution , to find the system at time in the spin state , with electrons and excess energy accumulated at the bath. One can readily write an equation of motion for this function, analogous to Eq. (27), to be Fourier transformed using four counting fields,

This quantity satisfies an equation of motion that is analogous to Eq. (29). It can be readily proved that the negative of the smallest eigenvalue of the corresponding matrix obeys the symmetry

which can be translated into the FT for the probability itself,

(47) | |||||

Here, . Enforcing energy and charge conservation, and , we recover Eq. (44).

### iii.4 Currents, and measures for vibrational cooling, heating, or instability

Currents. Analytical expressions for the charge and energy currents are obtained using the definition Eq. (33), utilizing Eqs. (LABEL:eq:mu) and (LABEL:eq:QE). These currents are defined positive when flowing to , and their closed forms are

(48) |

and

The TLS population is calculated in the steady-state limit,