First-passage time theory of activated rate chemical processes in electronic molecular junctions
aa r X i v : . [ c ond - m a t . m e s - h a ll ] J a n First-passage time theory of activated rate chemical processes in electronic molecularjunctions
Riley J. Preston
College of Science and Engineering, James Cook University, Townsville, QLD, 4811, Australia
Maxim F. Gelin
School of Sciences, Hangzhou Dianzi University, 310018 Hangzhou, China
Daniel S. Kosov
College of Science and Engineering, James Cook University, Townsville, QLD, 4811, Australia
Confined nanoscale spaces, electric fields and tunneling currents make the molecular electronicjunction an experimental device for the discovery of new, out-of-equilibrium chemical reactions.Reaction-rate theory for current-activated chemical reactions is developed by combining a Keldyshnonequilibrium Green’s functions treatment of electrons, Fokker-Planck description of the reactioncoordinate, and Kramers’ first-passage time calculations. The NEGF provide an adiabatic potentialas well as a diffusion coefficient and temperature with local dependence on the reaction coordinate.Van Kampen’s Fokker-Planck equation, which describes a Brownian particle moving in an externalpotential in an inhomogeneous medium with a position-dependent friction and diffusion coefficient,is used to obtain an analytic expression for the first-passage time. The theory is applied to severaltransport scenarios: a molecular junction with a single, reaction coordinate dependent molecularorbital, and a model diatomic molecular junction. We demonstrate the natural emergence of Lan-dauer’s blowtorch effect as a result of the interplay between the configuration dependent viscosityand diffusion coefficients. The resultant localized heating in conjunction with the bond-deformationdue to current-induced forces are shown to be the determining factors when considering chemicalreaction rates; each of which result from highly tunable parameters within the system.
I. INTRODUCTION
A molecular junction is a single molecule confined inthe nanoscale gap between two macroscopic, conductingleads. An applied voltage allows for the flow of electroniccurrent across the system through the valence states ofthe molecule. Large current densities and power dissipa-tion give way to strong current-induced forces and bond-selective heating which act to destabilize the molecularconfiguration within the junction, resulting in molecularconformational changes including telegraphic switching[1–4], along with providing the necessary energy for to-tal bond rupture[5–11]. This is obviously an undesirablefeature for promoting molecular electronics as a possibleavenue for moving beyond the traditional silicon semi-conductor technology into a regime of highly efficient andtailorable molecular-scale devices, and so a thorough the-oretical understanding is required for further progress.Nonetheless, it creates an exciting opportunity to exploreand produce new chemical reactions by providing a de-vice which traps a single molecule in a confined space ofa few nanometers where the electric field and current areapplied locally and selectively [12–14].The adequate and well established theories have beendeveloped for reaction rate calculations in gas and con-densed phases [15–22]. However, the development of sim-ilar theories for molecules in an electronic junction en-vironment is no simple task and as such, the scope oftheoretical work is still very limited. Three types of ap-proaches have been proposed to model current-induceddissociation. The first is based on the rate equation ap- proach where a single harmonic vibration is pumped be-yond the dissociation threshold limit [23, 24]. The sec-ond is a numerically exact scheme, which uses the hierar-chical quantum master equation method in conjunctionwith a discrete variable representation for the nucleardegrees of freedom to numerically study current-induceddissociation [25]. The third uses Keldysh nonequilibriumGreen’s functions to obtain a Fokker-Planck equation forthe reaction coordinate which is used to compute averageescape times and the accompanying reaction rates [5, 6].The further development of this approach is the subjectof this paper.Our approach utilizes the Born-Oppenheimer approx-imation, in which nuclei within the system are assumedto be slow-moving, classical particles interacting with asea of fast, quantum electrons. This separation of time-scales enables the use of a Langevin description for themotion of nuclei, in which the forces due to the quan-tum electrons are conceived through a frictional forcewhich acts to drive the nuclei into a motionless state,a fluctuating force which re-energizes the nuclei, and anadiabatic force which can deform the structure of the re-action potential. This enables the consideration of highlynon-trivial behaviour on the nuclear dynamics at thecost of a fully quantum description. Nevertheless, themethod has proven successful in a range of circumstances[1, 3, 4, 6, 11, 26–32].The method further lends itself to the study ofcurrent-induced chemical reaction rates [5, 6]. The useof Langevin dynamics to compute reaction rates wasfirst explored by Kramers in his seminal 1940 paper[33], in which the mean first-escape time of a parti-cle trapped in an arbitrary potential well subject toLangevin forces was approximated. The next significantstep was made in the 1990s, when Kramers’ theory wasextended to account for position-dependent friction andgeneralized Langevin equations describing finite-memory(non-Markovian) fluctuation-dissipation processes [34–37]. The effect of a velocity-dependent friction onKramers’ rates was investigated in Ref. [38]. However,these studies were limited to the regime of thermody-namic equilibrium. Beyond this regime, the fluctuation-dissipation theorem no longer holds, allowing for theemergence of localized heating effects in analogy to Lan-dauer’s proposed blowtorch effect [39, 40], in which spe-cific configurations of the reaction coordinate may expe-rience heightened temperatures, which may have a signif-icant affect on the evolution of the system. Such systemsare not limited to the realms of molecular electronics;the most common examples include numerous molecularmotors, ratchets, and heat engines [41–43] as well as var-ious confined nanosystems [44–48], notably of biologicalsignificance [49, 50]. Several explicit simulations of Lan-dauer’s blowtorch effects in double-well potentials havealso been performed recently [51–53].One of the aims of this paper is to further shed light onthis topic. A comprehensive understanding of the effectsof localized heating on the stability of molecular geome-tries is required to ensure the productive development ofspecific functionalities of molecular-scale electronic sys-tems. In this paper, we relax the requirement of ther-modynamic equilibrium, allowing for the self-consistentstudy of the mean first-passage times in model molecularelectronic junctions in both the underdamped and over-damped regimes. This is calculated through a Fokker-Planck equation, which arises due to our Langevin de-scription of the reaction coordinate within the junction.To a certain extent this work is the continuation oftwo papers of one of the authors (DSK) [5, 6], howeverRef.[5] considered the problem employing the fluctuation-dissipation theorem and Ref.[6] focused on the under-damped case only.The paper is structured as follows. In section II, wedemonstrate our calculations for the mean first-passagetime in the limiting regimes. This involves the calculationof the current-induced forces in the system, from which aFokker-Planck description then yields an equation for themean first-passage times. This is then applied to a simplemodel of the blowtorch effect in section III A. In sectionIII B, we calculate the reaction rates for a single-leveljunction model, in which current-induced forces are cal-culated self-consistently within the model. This is thenfurther applied to a model two-level molecule within thejunction in section III C. Atomic units are used in all cal-culations such that ~ = e = 1. Boltzmann’s constantis set to one k B = 1 in all derivations and calculationsmeaning that the temperature is measured in units ofenergy. II. THEORYA. Hamiltonian
We begin with the general Hamiltonian describing amolecular junction as given by H ( t ) = H M ( t ) + H L + H R + H LM + H MR . (1)The total system Hamiltonian is partitioned into the fol-lowing components; the molecular Hamiltonian H M ( t ),the left and right leads Hamiltonians H L and H R , and theleads-molecule coupling Hamiltonians H LM and H MR which describe the coupling between the electronic stateson the central molecule and the left and right leads, re-spectively.The molecular Hamiltonian takes the form: H M ( t ) = X ij h ij ( x ( t )) d † i d j + p m + U ( x ) , (2)where the operators d † i and d j denote the creation andannihilation operators for the molecular electronic stateswhose energies are given by the Hamiltonian matrix ele-ments h ij ( x ( t )). Note the explicit time dependence here,which arises as a result of the time evolution of the clas-sical reaction coordinate x . The last two terms in H M are not quantum mechanical operators, p / m is the ki-netic energy for the motion of the reaction coordinateand U ( x ) is the classical potential energy.The leads Hamiltonian is taken in the standard formof non-interacting electrons reservoirs: H L + H R = X kα ǫ kα d † kα d kα , (3)where the creation and annihilation operators are givenby d † kα and d kα , and the subscript kα denotes the op-erator acting on state k in the α lead which has energy ǫ kα .Finally, the system-lead coupling Hamiltonians H LM and H MR are given by: H LM + H MR = X kαi (cid:16) t kαi d † kα d i + h.c. (cid:17) . (4)The matrix elements t kαi (and their conjugates) describethe tunnelling amplitudes between lead states kα and themolecular orbitals i . B. Green’s Functions and Self-Energies
The foundational components from which we build ourmodel for the classical motion within a molecular junc-tion are the adiabatic Green’s functions. The requiredGreen’s function components (advanced, retarded, lesserand greater)[54] can be solved for from the Keldysh-Kadanoff-Baym equations, expressed in the Wigner spaceby[55] (cid:16) ω + i ∂ t − e i ∂ Gω ∂ ht h ( t ) (cid:17) G R/A ( t, ω ) = I + e − i ∂ Σ ω ∂ G t Σ R/A ( ω ) G R/A ( t, ω ) , (5)and (cid:16) ω + i ∂ t − e i ∂ Gω ∂ ht h ( t ) (cid:17) G > ( t, ω ) = e − i ∂ Σ ω ∂ Gt (cid:16) Σ R ( ω ) G > ( t, ω ) + Σ > ( ω ) G A ( t, ω ) (cid:17) , (6)where we have shown the retarded/advanced andlesser/greater equations collectively. Here, we have al-ready assumed that our molecule-leads coupling compo-nents t kαi are time-independent, which makes the self-energies only ω -dependent.In alignment with our previous work[3, 4, 32, 56–58],we assume that the classical motion along the reactioncoordinate within the system occurs over long time-scalesrelative to the characteristic electron tunnelling time.This provides us with the required small parameter tobe able to perturbatively solve (5) and (6) up to the firstorder in expansion of the exponents with derivatives. Re-sultantly, the adiabatic Green’s functions in the Wignerspace (Green’s functions which instantaneously follow thechanges in the reaction coordinate) are given by matricesof the form: G R/A ( x, ω ) = ( ωI − h ( x ) − Σ R/A ) − , (7) G > ( x, ω ) = G R ( x, ω )Σ > ( ω ) G A ( x, ω ) , (8)where I represents the identity operator in the molecu-lar space. The corresponding self-energy components aregiven by Σ Rα,ij = − i α,ij , (9)Σ Aα,ij = i α,ij , (10)Σ <α,ij ( ω ) = if α ( ω )Γ α,ij , (11)and Σ >α,ij ( ω ) = − i (1 − f α ( ω ))Γ α,ij , (12)where we have applied the wide-band approximation tothe leads, eliminating any ω dependence from Γ. As a result, the elements of the level broadening function aregiven by Γ α,ij = 2 πρ α t ∗ αi t αj , (13)where ρ α is the constant, energy-independent density ofstates for lead α and t iα is the tunneling amplitude t ikα which no longer depends on k .We additionally solve for the first non-adiabatic cor-rection to the Green’s functions, which account for thenon-zero velocity of the classical coordinate. The Green’sfunction corrections are given by [31] G R/A (1) = 12 i G
R/A h G R/A , ∂ t h i − G R/A , (14)and G < (1) = G R Σ < G A (1) + G R (1) Σ < G A + 12 i G R (cid:16) ∂ T hG R ∂ ω Σ < + G < ∂ T h + h.c. (cid:17) G A . (15) C. The Langevin Equation
Under the overarching assumption that the reactioncoordinate x along with its corresponding momentum p are classical variables in our approach due to the sep-aration of time-scales within the system, the equationof motion of the reaction coordinate can be expressedin the form of a quasi-classical Langevin equation, inwhich the classical motion is dictated by quantum me-chanical forces. Our Langevin equation is given by[1, 3, 4, 6, 30, 31] dpdt = − ∂ x U ( x ) + F ( x ) − ξ ( x ) m p + δf ( t ) , (16)in which we have the external classical potential U ( x ),an adiabatic force F ( x ) arising due to the occupation ofelectronic levels within the molecule, the frictional forceand its corresponding viscosity coefficient ξ ( x ), and thestochastic force δf ( t ). In our model, the adiabatic forcetakes the form F ( x ) = i π Z dω Tr n ∂ x hG < o . (17)Similarly, the viscosity coefficient is defined accordingto[6, 31] ξ ( x ) m p = − i π Z dω Tr n ∂ x hG < (1) o . (18)It is assumed that the fluctuating force can be de-scribed by a Gaussian process, defined by h δf ( t ) i = 0 , (19)and h δf ( t ) δf ( t ′ ) i = D ( x ) δ ( t − t ′ ) , (20)where D ( x ) is the quantum-mechanically calculated dif-fusion coefficient which determines the variance in thefluctuating force. In our model, D ( x ) is calculated ac-cording to[6, 31] D ( x ) = 12 π Z dω Tr n ∂ x hG > ∂ x hG < o . (21) D. Effective potential energy surface
The key quantity for chemical reaction rate calcula-tions is the effective potential energy surface. Here theeffective nonequilibrium potential energy surface is de-fined via the integration of the nonequilibrium force inthe Langevin equation U eff ( x ) = U ( x ) − Z xx dy F ( y ) . (22)We chose x in such a way that U eff ( x ) at the bottom ofthe potential well, x a , is zero, U eff ( x a ) = 0. This enforcesthat the energy E in the formulae below can vary in therange from 0 to ∞ . Suppose that the effective potentialhas the ”standard” Kramers problem form, in which casewe have a minimum at x = x a (reactant state) and anenergy barrier at x = x b , separating the reactant andproduct state.Let’s call U a = U eff ( x a ) = 0 (23)and U b = U eff ( x b ) . (24) E. Calculations of chemical reaction rates
1. Fokker-Planck equation with reaction-coordinatedependent viscosity and diffusion
The starting point for all derivations of the first-passage times is van Kampen’s Fokker-Planck equationfor the probability density function ρ ( x, p, t ) which de-scribes a Brownian particle of mass m , position x andmomentum p moving in the external potential U eff ( x ) inan inhomogeneous medium with position-dependent fric-tion ξ ( x ) and diffusion coefficient D ( x ) [59] ∂ t ρ ( x, p, t ) = (cid:16) − ∂ x pm + ∂ p U ′ eff ( x ) + n ∂ p pm ξ ( x ) + ∂ p D ( x ) o(cid:17) ρ ( x, p, t ) . (25)Van Kampen’s Fokker-Planck equation (25) is internallyconsistent (for example, it yields energy- and mass-flowtransport coefficients obeying Onsager’s reciprocity rela-tions [60]) and can unambiguously be derived from theLangevin equation (16) with Gaussian stochastic forceobeying (19) and (20). The Itˆo-Stratonovich dilemma[61] does not arise in the derivation of van Kampen’sFokker-Planck equation (25), because the viscosity ξ ( x )of Eq. (18) and the diffusion coefficient D ( x ) of Eq. (21)are momentum-independent.
2. Underdamped limit
In this section, we present our equations for the meanfirst-passage times from a reaction potential in an inho-mogeneous medium for both the underdamped and over-damped limiting regimes. Beginning with van Kampen’sFokker-Planck equation (25) for our reaction coordinate x and momentum p , we can obtain a one-dimensionalenergy-diffusion equation which is valid in the under-damped limit (see Appendix A for the derivation): ∂ t ρ ( E, t ) = ∂ E D ( E ) ρ unst ( E ) ∂ E ( ρ unst ( E )) − ρ ( E, t ) . (26)Here, we have introduced the energy-diffusion D ( E ) as D ( E ) = ν ( E )Ω( E ) , (27)defined as the ratio of the two energy-dependent func-tions, ν ( E ) = r m Z dxξ ( x ) T eff ( x ) p E − U eff ( x ) , (28)and Ω( E ) = r m Z dx p E − U eff ( x ) , (29)where we have defined the effective temperature T eff as T eff ( x ) = D ( x )2 ξ ( x ) . (30)In Eqs. (28) and (29) and in all other similar expres-sions, the integration domain corresponds to E > U eff ( x ).The stationary distribution is written as ρ unst ( E ) = Z − Ω( E ) exp ( − Z E dE ′ µ ( E ′ ) ν ( E ′ ) ) , (31)with normalization constant Z un = Z dE Ω( E ) exp ( − Z E dE ′ µ ( E ′ ) ν ( E ′ ) ) , (32)and µ ( E ) = r m Z dxξ ( x ) p E − U eff ( x ) . (33)The energy diffusion equation can be used to calculatethe mean first-passage time in the underdamped regimeas per τ = Z U b U a dE ′ D ( E ′ ) ρ unst ( E ′ ) Z E ′ dE ′′ ρ unst ( E ′′ ) . (34)
3. Overdamped limit
In the overdamped limit, van Kampen’s Fokker-Planckequation (25) reduces to the diffusion equation [59, 62]which can be cast into a form similar to Eq. (26): ∂ t ρ ( x, t ) = ∂ x T eff ( x ) ξ ( x ) ρ ovst ( x ) ∂ x ( ρ ovst ( x )) − ρ ( x, t ) . (35)The stationary distribution is given by ρ ovst ( x ) = 1 Z ov T eff ( x ) exp (cid:26) − Z x dx ′ ∂ x U ( x ′ ) − F ( x ′ ) T eff ( x ′ ) (cid:27) , (36)where Z ov = Z + ∞−∞ dx T eff ( x ) exp (cid:26) − Z x dx ′ ∂ x U ( x ′ ) − F ( x ′ ) T eff ( x ′ ) (cid:27) . (37)The first-passage time in the overdamped regime canbe thus written as τ = Z x b x a dx ′ ξ ( x ′ ) T eff ( x ′ ) ρ ovst ( x ′ ) Z x ′ −∞ dx ′′ ρ ovst ( x ′′ ) . (38)To summarize the derivations in section II: the mainresult is the combined use of the adiabatic force (17),the position-dependent viscosity (18) and the position-dependent diffusion coefficient (21) - all obtained fromnonequilibrium Green’s function calculations - with ex-act (in underdamped and overdamped limits) expressions(34,38) for the mean first-passage time τ . The rates forcorresponding chemical reactions can obtained as inverseof the first-passage time 1 /τ . III. RESULTSA. The blowtorch effect
In this section we investigate Landauer’s proposedblowtorch effect[40], in which a non-equilibrium systemallows for coordinate-dependent variations to the dissi-pative forces acting on a particle, which then has an ef-fect on the properties of the steady-state distribution.Landauer’s blowtorch effect plays a critical role in chem-ical reactions in molecular electronic junctions, thereforewe first discuss its general features which will be rele-vant for our subsequent discussion. For consistency withKramers’ seminal paper[33] on chemical reaction rates,our analysis is formulated in terms of the mean first es-cape time from the left minimum of a bistable potential,which we model according to a quartic of the form U ( x ) = − ax + bx , (39)where x is our reaction coordinate. The constants a and b are adjustable parameters which determine thewidth and depth of the minimum. In all tests in this sec-tion, we set a = 0 .
04 and b = 0 . U b = 0 . ξ and diffusion coefficients D are set to a constant value over the range of x , yieldinga constant temperature as determined by the fluctuation-dissipation theorem.In order to introduce an inhomogeneity into the tem-perature, a Gaussian spike is applied to the diffusion co-efficient locally at position x , D ( x ) = D + D peak ( x ) , (40)where D peak ( x ) = D m e − ( x − x σ , (41)with adjustable width σ and magnitude D m parameters.The effective temperature profile is given by (30), theeffective temperature at the peak is then given by T max = T (1 + D m D ) , (42)where T = D ξ . (43)This represents Laundauer’s so-called ”blowtorch”which heats a small segment of the reaction coordinate,as shown diagramatically in Fig.1.FIG. 1: An adjustable temperature spike is introducedwhich heats a chosen part of our reaction coordinatepotential.Here, the intention is to study the effect of shifting thetemperature spike along the reaction coordinate on themean first-passage time τ . We analyse the overdampedand underdamped regimes separately for the same pa-rameters except for the mass m of a Brownian particle,which is chosen to satisfy the desired regime.In Fig.2, we observe the behaviour of the mean first-passage time as the position of the temperature peak isshifted along the reaction coordinate (shown in blue),while the reaction potential is shown as a reference inorange. In the underdamped regime, we observe τ tobe minimized when the heating is applied to the bottomof the potential. This enables the molecule to heat upquickly at low energies, and repeatedly attain more en-ergy as it passes through this region in a near-harmonicmanner.The overdamped regime differs, in that τ is minimizedwhen the heating is applied approximately halfway upthe potential, around the point of steepest ascent. Inthe overdamped regime, the escaping particle will veryquickly equilibrate to any given temperature fluctuationto which it is exposed. As such, the heated region causesa flattening of the probability distribution in that re-gion, nullifying the dependence of the distribution onthe reaction coordinate. This causes an effective reduc-tion to the height of the energy barrier U b as elucidatedby Landauer[39, 40]; a phenomenon which is maximizedwhen the heating is applied in the region of steepest as-cent. We note the counter-intuitive observation that ifthe heating is applied to the bottom of the potential inthe overdamped case, this causes only a small reductionto τ . This is because the particle will quickly lose theobtained energy as it returns to the cooler regions whenit attempts to escape.It is insightful for us to also study the effect of thestrength of interaction of a Brownian particle with theenvironment, while maintaining a homogeneous temper- -0.5 0 0.5 1 1.51234567 10 (a) -0.5 0 0.5 1 1.50.511.52 10 (b) FIG. 2: The calculated mean first-passage time τ as afunction of the temperature peak’s position along thereaction coordinate, for the (a) overdamped and (b)underdamped regime ( m = 1000a.u). The black linedenotes τ in the absence of an applied blowtorch.Parameters: D = 0 . ξ = 1, σ = 0 . x will be counteracted by acorresponding changes in the viscosity coefficient at thesame x , enforcing a homogeneous temperature as per thefluctuation-dissipation theorem. Here, we perform a sim-ilar analysis as above, such that we have a moveable peakof increased interaction (simultaneously locally increaseddiffusion and viscosity) while the temperature is homo-geneous. The results of this are displayed in Fig.3.In the underdamped regime, we observe that thelargest reduction to τ occurs when the interaction peakis placed at the minimum. The decrease in τ agrees withthe homogeneous-case solution, with the distinction thatreaction coordinates at higher energies in the potentialhave diminishing contributions to the decreasing τ . In -0.5 0 0.5 1 1.5678910 10 (a) -0.5 0 0.5 1 1.51.851.91.9522.052.12.15 10 (b) FIG. 3: The calculated mean first-passage time as afunction of the interaction peak’s position along thereaction coordinate, for the (a) overdamped and (b)underdamped regime ( m e = 1000a.u). The black linedenotes τ in the absence of an applied blowtorch.Parameters: D = 0 . ξ = 1, σ = 0 . τ as also predicted in thehomogeneous case. However, we observe that this is dom-inated by the increased interaction strength near to themaximum of the potential, while changing the interac-tion strength in the rest of the potential has negligibleeffect. This demonstrates that τ has little regard for theinteraction strength in the overdamped regime, except inthe region approaching the maximum.This general analysis arms us with the required phys-ical intution before proceeding to the next section, inwhich we first observe how a Landauer blowtorch emergesnaturally from a simple molecular junction model, thendemonstrate the effect on hypothetical chemical reactionrates. B. Application to a molecule with a singlecurrent-carrying molecular orbital
In this section, we analyze the calculated mean first-passage time τ for our model of a molecular junction.Contrary to section III A, the viscosity and diffusion co-efficients will be computed using nonequilibrium Green’sfunctions according to eqs. (18) and (21). We considerthe case of a single electronic level coupled to the leftand right leads under some applied bias voltage. Themolecular Hamiltonian is simplified to H M = ( h + λx ) d † d, (44)where subscript d † and d † denote the creation and an-nihilation operators for an electron on the molecular or-bital, while all previous matrix quantities are simplifiedto scalars in this regime. Here, the dependency of H M onthe reaction coordinate acts to shift the electronic level,as scaled by the tuneable parameter λ . The left and rightlead are each at room temperature (0.00095a.u) and aresymmetrically coupled to the central electronic state suchthat our level-width function is given by Γ L = Γ R = 0 . x a : U ( x ) = − a ( x − x a ) + b ( x − x a ) . (45)This means that when x a = 0, the potential mini-mum (ignoring the effects of F ) occurs at x = 0, whilea positive x a shifts the input potential to the right.Any bias voltage is applied symmetrically, such that µ L = − µ R = V /
2, where µ L and µ R are the chemicalpotentials of the left and right leads.We study the effect of applying a gate voltage to thesystem, as modelled by a shift in the h value. This al-lows for a degree of controllability of the reaction ratesfor a given system. Figs.4a and 4b demonstrate the re-sultant viscosity coefficient and temperature respectively,as a function of the reaction coordinate. Application of agate voltage acts to shift the curve along the reaction co-ordinate. This analysis is performed for a non-zero biasvoltage such that the temperature is now inhomogeneousin addition to the viscosity. In the underdamped regimeshown in Fig.4d, τ is minimized when the viscosity andtemperature peaks are shifted near to the minimum of thereaction potential (note however, that the minimum in τ does not occur exactly when the peaks are shifted to theminimum of the potential due to the slight asymmetry -5 0 500.010.020.03 (a) -5 0 555.56 10 -3 (b) -4 -2 0 2 4 600.511.522.5 10 (c) -4 -2 0 2 40246810 10 (d) FIG. 4: The effect of an applied gate voltage to (a) the viscosity coefficient and (b) the effective temperature. Themean first-passage time τ in the (c) overdamped and (d) underdamped regime is plotted against the peak coordinateof the viscosity and temperature (as determined by the applied gate voltage) for different λ . The coordinates of theminimum and maximum of the reaction potential are denoted by the vertical black lines in (c) and (d).of the reaction potential). In contrast, the overdampedregime displays highly non-trivial behaviour, arising asa result of the interplay between the strength of the vis-cosity and the temperature. In our analysis of the over-damped regime in the previous section, we noted thatthe dependence of τ on the temperature is dominated bythe region of steepest ascent up towards the maximum.Here, we again observe this behaviour as the large peak inFig.4c corresponds to when the dip in the temperatureoccurs in this region (when the temperature peak hasbeen shifted to the right). A corresponding but smallerpeak also occurs due to a shift to the left in the tempera-ture such that the low temperature aligns with the steepregion of the potential. The difference in peak sizes arisesas a result of the inhomogeneous viscosity, which per theprevious section, we know is important in the region nearthe maximum of the reaction potential. The large peakin τ occurs when the temperature is low in the steep re-gion, while the viscosity is high towards the maximum.The small peak has a low viscosity near the maximum,explaining its comparatively smaller magnitude. C. Model of a two-level molecule
In this section, we expand the model to consider atwo-level system. The second molecular orbital is not asimple addition of an extra level here, since the Green’sfunctions, self-energies and molecular Hamiltonian be-come 2 × H +2 molecule[63]. As such, the molecular Hamiltonian nowreads H M ( t ) = X ij h ij ( q ( t )) d † i d j , (46)where d † i and d j are now in the molecular orbital basis.The electronic Hamiltonian elements are represented inthe form of a 2 × h = (cid:18) H b ( q ) 00 H a ( q ) (cid:19) , (47)where we use H b ( q ) and H a ( q ) to denote the bondingand anti-bonding molecular orbitals, respectively, while q is the bond-length. The values for H b ( q ) and H a ( q ) arecalculated according to molecular orbital theory[63] andare given by H b ( q ) = H AA + H AB S AB , (48)and H a ( q ) = H AA − H AB − S AB , (49)where H AA and H AB are the Hamiltonian elements inthe atomic basis and S AB is the overlap integral betweenatomic 1s Slater orbitals. The constituent componentsare then given by H AA = −
12 + e − q (cid:16) q (cid:17) − q , (50) H AB = − S AB ( q )2 − e − q (1 + q ) , (51)and S AB = e − q (1 + q + q / . (52)In the interest of simplicity, each of the molecular or-bitals is symmetrically coupled to the left and right leads;as controlled by the parameter Γ which now takes theform of a matrix as perΓ α = (cid:18) Γ α,bb Γ α,ba Γ α,ab Γ α,aa (cid:19) (53)for the α lead, where the off-diagonal components can bedefined according to,Γ α,ba = Γ α,ab = p Γ α,aa Γ α,bb . (54)In each test, we have µ L = − . µ R = µ L − V , whilethe lead temperatures are again set to room temperature.The external potential now represents the classical nu-clear repulsion, which in atomic units is given by U ( q ) = 1 q . (55) Inclusion of the electronic forces allows us to generatemodified electronic potentials for varied parameters inorder to assess the molecular stability. Examples of thesepotentials are shown in Fig.5; where in (a) an appliedbias voltage is shown to decrease the energy required forbond rupture, while (b) shows the effect of the additionalelectronic level which when occupied, acts to increase thebond stability.Along with the shape of the effective potential, thebond stability is also determined by the electronic vis-cosity and effective temperature, which are demonstratedin Fig.6. In the viscosity coefficient, each curve showsa peak at small q , which approximately corresponds towhen H a crosses the fermi-level of the left lead. Likewise,the peaks at large q are a result of H b crossing the fermi-level of the right, then left, leads (these split peaks mergetogether when V = 0). The inset plot demonstrates theeffect of allowing an additional transport channel throughthe excited state which not only introduces the peak atsmall q , but also increases the magnitude of the viscousforces overall. The effective temperature is equal to theleads temperature for V = 0, while non-zero voltagesyield a complex array of localized heating and cooling ef-fects, which arise as the energy levels shift in and out ofthe resonance region as the bond-length is increased.These competing effects culminate in our calculationof the mean first-passage time, which is demonstrated inFig.7 as a function of the bias voltage, for different cou-pling values to the excited electronic state. In both limit-ing regimes, an increase to the bias voltage acts to desta-bilize the bond and decrease τ , both due to the increasedeffective temperatures and the weakening of the bonddue to the current-induced forces. In the overdampedregime, allowing the leads to be coupled to an additionallevel in the central region has a stabilizing effect for allvoltages tested, increasing the average amount of timefor bond rupture. The underdamped case shows similarbehavior for very low voltages; however, at higher volt-ages the availability of the additional transport channelthrough the excited state increases the current-inducedforces such that the energy required for a bond-ruptureis found more easily, decreasing τ . IV. CONCLUSIONS
In this paper, we have demonstrated that the rates ofchemical reactions for molecules in electronic junctionsdepend on three crucial ingredients; the potential energysurface which defines the energy required for a configu-ration change or bond rupture, the rate of the energy re-moval from vibrational to electronic degrees of freedomgiven by the electronic viscosity coefficient, and lastly,the effective temperature dynamically established in themolecule. While the magnitude of these quantities is ofhigh importance, the local distribution of the viscosityand effective temperature along the potential energy sur-face (Landauer’s blowtorch effect) also proves to be crit-0 (a) (b)
FIG. 5: The adiabatic potential as a function of the bond-length presented for (a) varying bias voltages and (b)varying the magnitude of leads coupling to H a . Parameters: V = 0, Γ aa / Γ bb = 1, unless otherwise specified.Γ bb = 0 .
03 in all calculations. (a) -3 (b) FIG. 6: The effect of varying the bias voltage is shown for the (a) viscosity coefficient and (b) the effectivetemperature, as a function of the bond length. Insets: Shows the same quantity at V = 0 .
02 for Γ aa / Γ bb = 0(dashed) and Γ aa / Γ bb = 0 . aa / Γ bb = 1 in the main plots. Γ bb = 0 .
03 in all calculations. (a) (b) FIG. 7: The mean first-passage time τ as a function of the bias voltage, varying the coupling to H a in the (a)overdamped and (b) underdamped cases.1ical.The addition of localized heating and cooling effects asa result of inhomogeneity with respect to the molecularconfiguration has been shown to induce significant varia-tions in the mean first-passage time, as calculated accord-ing to a Fokker-Planck description obtained by separat-ing slow (reaction coordinate) and fast (electronic) time-scales in the Keldysh-Kadanoff-Baym equations for thenonequilibrium Green’s functions. This has been demon-strated for a single-level molecular junction model, aswell as a two-level model inspired by H +2 molecular or-bitals with the bond length considered as the reactioncoordinate. This interplay between the amount of en-ergy required for bond rupture and the energy supplieddue to tunneling electrons has been shown to be stronglydependent on the choice of experimentally tuneable pa-rameters for the system. This enables the possibility of ahigh degree of controllability for molecular junction sys-tems, with promises of controlled initiation of chemicalreactions or conversely, enforcing the stability of specificconfigurations within the system. DATA AVAILABILITY
The data that supports the findings of this study areavailable within the article.
Appendix A: Derivation of the energy diffusionequation in the underdamped limit
Following Zwanzig [15], we can introduce the projec-tion operator b O = Ω − ( E ) Z dxdp δ ( E − H ) (A1)where H = p m + U eff ( x ) (A2)is the Hamiltonian of the bath-free Brownian particleand the microcanonical partition function is defined asΩ( E ) = Z dxdp δ ( E − H ) (A3)(equivalently, Ω( E ) is determined via Eq. (29)). Notethat b O has an extra factor of Ω − ( E ) in comparison withthe operator introduced by Zwanzig. With our defini-tion, b O = b O . Applying the projection operator of Eq.(A1) to the Fokker-Planck equation (25), exploiting theidentity b O (cid:16) − ∂ x pm + ∂ p U ′ eff ( x ) (cid:17) = 0and taking into account that ξ ( x ) ≪ ∂ t b Oρ ( x, p, t ) = b Oζ ( x ) n ∂ p pm + ∂ p T ( x ) o b Oρ ( x, p, t ) . (A4) Changing the differentiation variables ( ∂ p = m − ∂ E p )and using the identity ∂ p b O = m − ∂ p p∂ E b O , we obtainthe following equation for ρ ( E, t ) ≡ Ω( E ) b Oρ ( t ) (cf. [15,38]): ∂ t ρ ( E, t ) = ∂ E { µ ( E ) + ν ( E ) ∂ E } Ω − ( E ) ρ ( E, t ) . (A5)Here µ ( E ) and ν ( E ) are defined through Eqs. (33)and (28), and ρ ( E, t ) is normalized according to R dEρ ( E, t ) = 1.2 [1] G. Weick, F. Pistolesi, E. Mariani, and F. von Oppen,Phys. Rev. B , 121409 (2010).[2] J.-T. L¨u, M. Brandbyge, P. Hedeg˚ard, T. N. Todorov,and D. Dundas, Phys. Rev. B , 245444 (2012).[3] R. J. Preston, V. F. Kershaw, and D. S. Kosov,Phys. Rev. B , 155415 (2020).[4] V. F. Kershaw and D. S. Kosov,The Journal of Chemical Physics , 154101 (2020),https://doi.org/10.1063/5.0023275.[5] A. A. Dzhioev and D. S. Kosov,The Journal of Chemical Physics , 074701 (2011),https://doi.org/10.1063/1.3626521.[6] A. A. Dzhioev, D. S. Kosov, and F. von Oppen,The Journal of Chemical Physics , 134103 (2013),https://doi.org/10.1063/1.4797495.[7] A. Erpenbeck, C. Schinabeck, U. Peskin, and M. Thoss,Phys. Rev. B , 235452 (2018).[8] R. H¨artle and M. Thoss,Phys. Rev. B , 125419 (2011).[9] H. Li, T. A. Su, V. Zhang, M. L. Steiger-wald, C. Nuckolls, and L. Venkataraman,Journal of the American Chemical Society , 5028 (2015),pMID: 25675085, https://doi.org/10.1021/ja512523r.[10] H. Li, N. T. Kim, T. A. Su, M. L. Steigerwald, C. Nuck-olls, P. Darancet, J. L. Leighton, and L. Venkataraman,Journal of the American Chemical Society , 16159 (2016),pMID: 27960303, https://doi.org/10.1021/jacs.6b10700.[11] J.-T. L¨u, M. Brandbyge, andP. Hedeg˚ard, Nano Letters , 1657 (2010),https://doi.org/10.1021/nl904233u.[12] K. Huang, L. Leung, T. Lim, Z. Ning, and J. C. Polanyi,Journal of the American Chemical Society , 6220 (2013).[13] A. C. Aragon´es, N. L. Haworth, N. Darwish, S. Ciampi,N. J. Bloomfield, G. G. Wallace, I. Diez-Perez, and M. L.Coote, Nature , 88 (2016).[14] B. Borca, T. Michnowicz, R. P´etuya, M. Pristl,V. Schendel, I. Pentegov, U. Kraft, H. Klauk,P. Wahl, R. Gutzler, A. Arnau, U. Schlickum, andK. Kern, ACS Nano , 4703 (2017), pMID: 28437066,https://doi.org/10.1021/acsnano.7b00612.[15] R. Zwanzig, Nonequilibrium statistical mechanics (Ox-ford University Press, Oxford, 2001).[16] W. T. Coffey, Y. P. Kalmykov, and J. T. Waldron,
TheLangevin Equation (World Scientific, Singapore, 2004).[17] P. H¨anggi, P. Talkner, and M. Borcovec,Reviews of Modern Physics , 251 (1990),https://doi.org/10.1063/1.527963.[18] V. I. Melnikov, Phys. Rep. , 1 (1991),https://doi.org/10.1063/1.527963.[19] G. A. Voth and R. M. Hochstrasser,J. Phys. Chem. , 13034 (1996),https://doi.org/10.1063/1.527963.[20] W. H. Miller, J. Phys. Chem. A , 793 (1998),https://doi.org/10.1063/1.527963.[21] E. Pollak and P. Talkner, Chaos , 026116 (2005),https://doi.org/10.1063/1.527963.[22] B. Sch¨uller, A. Meistrenko, H. Hees, Z. Xu,and C. Greiner, Annals of Physics , 168045 (2020),https://doi.org/10.1063/1.527963.[23] J. Koch, M. Semmelhack, F. von Oppen, and A. Nitzan,Phys. Rev. B , 155306 (2006). [24] D. Brisker and U. Peskin,The Journal of Chemical Physics , 244709 (2008),https://doi.org/10.1063/1.3021288.[25] A. Erpenbeck, Y. Ke, U. Peskin, and M. Thoss,Phys. Rev. B , 195421 (2020).[26] W. Dou and J. E. Subotnik,Phys. Rev. B , 064303 (2018).[27] W. Dou and J. E. Subotnik,Phys. Rev. B , 104305 (2017).[28] W. Dou, G. Miao, and J. E. Subotnik,Phys. Rev. Lett. , 046001 (2017).[29] N. Bode, S. V. Kusminskiy, R. Egger, and F. von Oppen,Phys. Rev. Lett. , 036804 (2011).[30] F. Pistolesi, Y. M. Blanter, and I. Martin,Phys. Rev. B , 085127 (2008).[31] N. Bode, S. V. Kusminskiy, R. Egger, and F. von Oppen,Beilstein J. Nanotechnol , 144 (2012).[32] R. J. Preston, T. D. Honeychurch, and D. S. Kosov,The Journal of Chemical Physics , 121102 (2020),https://doi.org/10.1063/5.0019178.[33] H. Kramers, Physica , 284 (1940).[34] E. Pollak and A. M. Berezhkovskii,J. Chem. Phys. , 1344 (1993),https://doi.org/10.1063/1.527963.[35] G. R. Haynes, E. Pollak, and G. A.Voth, J. Chem. Phys. , 7811 (1994),https://doi.org/10.1063/1.527963.[36] G. R. Haynes and G. A. Voth,J. Chem. Phys. , 10176 (1995),https://doi.org/10.1063/1.527963.[37] E. Neria and M. Karplus,J. Chem. Phys. , 10812 (1996),https://doi.org/10.1063/1.527963.[38] M. F. Gelin and D. K. Kosov,J. Chem. Phys. , 244501 (2007),https://doi.org/10.1063/1.527963.[39] R. Landauer, Phys. Rev. A , 636 (1975).[40] R. Landauer, Physica A: Statistical Mechanics and its Applications , 551 (1993).[41] P. Reimann, Phys. Rep. , 57 (2002),https://doi.org/10.1063/1.527963.[42] P. H¨anggi and F. Marchesoni,Rev. Mod. Phys. , 387 (2009),https://doi.org/10.1063/1.527963.[43] S. Erbas-Cakmak, D. A. Leigh, C. T. McTernan,and A. L. Nussbaumer, Chem. Rev. , 10081 (2015),https://doi.org/10.1063/1.527963.[44] M. Das and D. S. Ray, Phys. Rev. E , 052133 (2015),https://doi.org/10.1063/1.527963.[45] J. Devine and M. W. Jack,Phys. Rev. E , 062130 (2017),https://doi.org/10.1063/1.527963.[46] G. Franzese, I. Latella, andJ. M. Rubi, Entropy , 507 (2017),https://doi.org/10.1063/1.527963.[47] V. Holubec, A. Ryabov, M. H. Yaghoubi,M. Varga, A. Khodaee, M. E. Foulaad-vand, and P. Chvosta, Entropy , 119 (2017),https://doi.org/10.1063/1.527963.[48] A. Arango-Restrepo and J. M. Rubi,J. Chem. Phys. , 034108 (2020),https://doi.org/10.1063/1.527963. [49] S. Basak, S. Sengupta, and K. Chat-topadhyay, Bioph. Rev. , 851 (2019),https://doi.org/10.1063/1.527963.[50] J. M. Rubi, EPL , 10001 (2019),https://doi.org/10.1063/1.527963.[51] M. Bekele, S. Rajesh, G. Ananthakrishna,and N. Kumar, Phys. Rev. E , 143 (1999),https://doi.org/10.1063/1.527963.[52] M. Das, D. Das, D. Barik, andD. S. Ray, Phys. Rev. E , 052102 (2015),https://doi.org/10.1063/1.527963.[53] M. Das and H. Kantz, Phys. Rev. E , 032108 (2019),https://doi.org/10.1063/1.527963.[54] G. Stefanucci and R. van Leeuwen, NonequilibriumMany-Body Theory of Quantum Systems: A Modern In-troduction (Cambridge University Press, 2013).[55] H. Haug and A. Jauho,