{'fraction_non_alphanumeric': 0.027471146154725175, 'fraction_numerical': 0.008958976117015198, 'mean_word_length': 5.107761027359017, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Non-Adiabatic Effects of Nuclear Motion in Quantum Transport of Electrons: A Self-Consistent Keldysh-Langevin Study

The molecular junction geometry is modelled in terms of nuclear degrees of freedom that are embedded in a stochastic quantum environment of non-equilibrium electrons. Time-evolution of the molecular geometry is governed via a mean force, a frictional force and a stochastic force -- forces arising from many electrons tunnelling across the junction. Conversely, the current-driven nuclear dynamics feed back to the electronic current, which can be captured according extended formula for the current that has an explicit dependency on classical nuclear velocities and accelerations. Current-induced nuclear forces and the non-adiabatic electric currents are computed using non-equilibrium Green's functions via a time-scale separation solution of Keldysh-Kadanoff-Baym equations in Wigner space. Applying the theory to molecular junctions demonstrated that non-adiabatic corrections play an important role when nuclear motion is considered non-equilibrium and, in particular, showed that non-equilibrium and equilibrium descriptions of nuclear motion produce significantly different conductivities. Finally, we observe that non-equilibrium descriptions of nuclear motion can give rise to the Landauer blowtorch effect via the emergence of multi-minima potential energy surfaces in conjunction with non-uniform temperature profiles. The Landauer blowtorch effect and its impact on the current characteristics, waiting times and the current noise are explored for an effective adiabatic potential morphs between a single, double and triple potential as a function of voltage.


I. INTRODUCTION
The transport of electrons through molecular-sized electronic systems leads to an abundance of new phenomena that are of high physical significance 1,2 . An important example of which is various phenomena related to the interaction between current and molecular vibrational degrees of freedom: a process that results in altered current characteristics of nanoscale electronic devices and, conversely, changed molecular geometries relative to their equilibrium configuration. Electrons tunnelling across the electronic device can exchange energy with the molecule as they travel across the junction; molecular vibrational modes can be activated [2][3][4][5][6][7][8][9] , often leading to extreme conformational changes in the molecular geometry, and possibly bond rupture within the molecule or at the molecule-lead interface [10][11][12][13] . With this, the changing nuclear degrees of freedom of the molecule may influence the transport properties of the device which, in turn, can lead to interesting transport phenomena such as negative differential resistance, Frank-Condon blockade and non-renewal statistics, amongst others 4,14-21 .
Experiment and theory have made leaps and bounds in obtaining a fundamental understanding of the interplay between electronic and nuclear degrees of freedom and its impact on electronic current. In terms of quantum-mechanical approaches, a common procedure is the use of non-equilibrium Green's functions (NEGF) to construct a various kinds of perturbative expansions about the strength of electron-vibrational coupling [22][23][24]24,25,25,26,26,27,[27][28][29][30][30][31][32][33][34][35][36][37][38][39][40][41][42][43][44] . Strong electronphonon interactions are usually treated according to master equation methods 17,20,21,[45][46][47][48][49][50][51][52][53][54][55] ; however, such formulations come at the cost of assuming that the molecule-metal coupling is small and, consequently, put it out of range of a large number of interesting molecular systems. All these methods, while providing tremendous insights into the workings of current-induced phenomena, suffer from one critical flaw: they restrict nuclear dynamics to vibrations about the fixed equilibrium geometry. According to these models, molecules in nanoscale molecular devices cannot undergo large-scale structural and/or conformational changes: the result is the neglect of important physical processes, which can have profound effects on the electronic transport properties of the considered devices, especially at the high voltage and strong electron-vibration coupling regimes.
An alternative approach to the problem of vibrationally-induced current is to re-envision the electronic coupling with nuclear dynamics according to that of 'electronic friction' acting on classical nuclei 13,[56][57][58][59][60][61][62][63][64][64][65][66][67][68][69][70][71] . According to this view, mechanical degrees of freedom of the molecule are modelled in terms of classical particles that are embedded in a stochastic electronic environment: time-evolution of the nuclear degrees of freedom is governed via a mean force, a frictional force and a stochastic force, forces arising from many electrons tunnelling across the junction for a given nuclear vibration. Conversely, the molecular junction geometry has impacts on the electronic current, which can be captured according to extended expression for the current that have explicit dependencies on classical nuclear velocities 72 and higher-order terms [67][68][69] . Such a formulation holds advantages when compared to some of the former theories mentioned before: it allows for easy descriptions of molecular geometries that can undergo radical deformations; furthermore, such a model provides a gateway into the description of non-adiabatic motion effects in the transport properties of molecular junctions as well as current-induced chemical reactions.
Of course, these benefits come at the cost of approximating quantum objects as classical, along with the restriction that nuclear motion must be considered 'slow' relative to the 'fast' tunnelling electrons.
This approach, which we will refer to henceforth as a Keldysh-Langevin description, has received renewed interest in recent years, particularly in the context of nonequilibrium electronic transport 13,[56][57][58][59][60][61][62]71 . In our previous works [67][68][69] , we developed extensions to the Meir-Wingreen formula 73 to consider terms that are second order in the central-time derivatives and demonstrated the existence of terms that are appreciable and nonvanishing when statically averaged over nuclear motion. In this paper, we extend this model by replacing the Boltzmann equilibrium description by a Langevin nonequilibrium description to explore how non-equilibrium nuclear dynamics affects the current characteristics of the device. The theory is based on the assumption that the characteristic time-scales associated with conformational changes of the molecule are slow relative to the fast tunnelling electrons, from which we construct a perturbative transport theory. While perturbative, our theory need not assume that molecular vibrations are small and harmonic about the equilibrium geometry, nor does it assume that the coupling between the electronic and nuclear degrees of freedom be considered small. The result of the theory is an extended formula for the electric current, along with a Langevin equation that selfconsistently simulates the current driven nuclear dynamics and current. This paper in many ways serves as an extension to our previous work [67][68][69] where nuclear degrees of freedom are no longer assumed to be in equilibrium but rather change dynamically due to stochastic forces exerted by nonequilibrium electrons. Section (II) of the paper illustrates the theoretical model: the Hamiltonian, self-energy and Green's functions are defined with a perturbative solution being implemented for the Keldysh-Kadanoff-Baym equations; and a set of Langevin equations and extended Meir-Wingreen formulas are derived. Section (III) then applies the proposed theory of electron transport to several problems in order to discuss new phenomena arising out of Langevin descriptions of nuclear dynamics. Note that this paper will work with the atomic unit convention of e = = m e = 1.

A. Hamiltonian
We start with the general set-up of molecular electronics for a nanoscale quantum system coupled to two macroscopic leads, a set-up that is governed by the Hamiltonian: Here the HamiltonianĤ(t) is the total system Hamiltonian and is comprised of a molecular HamiltonianĤ M (t), the left and right leads HamiltoniansĤ L andĤ R , along with the HamiltoniansĤ LM andĤ MR that describe the coupling between the molecule and the left and right leads, respectively. Note that we have made the time parameter t explicit to show which Hamiltonians carry explicit time-dependence.
Our study will consider a molecular Hamiltonian that has explicit time-dependence through a multidimensional vector x(t) which describes the nuclear degrees of freedom associated with the molecular geometry. The HamiltonianĤ M (t) takes the form: As can be seen above, the second-quantisation operatorsd † i andd j represent the creation and annihilation operators for the quantum single-particle states i and j, along with their respective Hamiltonian matrix elements h ij (x(t)).
The Hamiltonians for the leads are taken in the standard way as macroscopic reservoirs of non-interacting electrons: The creation and annihilation operatorsd † kα andd kα (for α ∈ {L, R}) create and annihilate an electron in the leads in the single-particle state k in the α th lead.
Finally, we have the system-lead coupling Hamiltoni-ansĤ LM andĤ MR that take the form: The matrix elements h kα,i (and their conjugates) describe the tunnelling amplitudes between lead states kα and the molecular orbitals i.

B. Green's Functions and Self-Energies
Given our time-independent coupling constants h kα,i and h i,kα , the self-energy depends on relative time only and therefore becomes time-independent in energy space. The retarded and advanced components take the form: and The level-width functions are given by where the level-shift functions ∆ αij (ω) can be computed from Γ αij (ω) via Kramers-Kronig relation. The lesser and greater self-energies are computed as and All our physical quantities including the non-adiabatic corrections will be expressed in terms of readily-available adiabatic Green's functions. All components of the Green's functions depend instantaneously on the molecular geometry and are defined though standard relations. The Green's function components are given by: and Above we have the identity matrix I in molecular orbital space.

C. The Langevin Equation
The nuclear degrees of freedom are considered as classical variables within our approach. If, additionally, we assume that nuclear dynamics are slow relative to the tunnelling electrons, we can obtain a Langevin equation of motion for each component µ 13,56,57,[60][61][62]71 , where we have used the notation ∂ µ ≡ ∂ ∂xµ . The Langevin equations above are comprised of a classical and external classical potential U (x), an adiabatic drag force F µ , a frictional force and its viscosity tensor ξ µν , and a stochastic force δf (t). The adiabatic drag force F µ has the form with the matrix Λ µ being the derivative of single-particle molecular Hamiltonian matrix: The viscosity tensor which depends on the molecular junction geometry and is given by 61 The stochastic force δf µ (t) is modelled as a Markovian Gaussian variable with mean δf µ (t) = 0 (18) and delta-function variance Above we have the diffusion coefficient D µν ; one can show that D µν can be expressed in terms of lesser and greater Green's functions 61 :

D. Electric Current
A derivation of a formula for electric current proceeds routinely by transforming the equations of motion to the Wigner space and expanding the exponential operator to consider only those terms that are second-order in central-time derivatives, an expansion that is performed in the same assumption as in the derivation of the Langevin equation for nuclear degrees of freedom that are slow relative to the tunnelling electrons. [67][68][69] The expression for the electric current can be cast into the form of the Landauer formula 67 where ... x means the statistical averaging over the Langevin trajectory. The subscript "pert" indicates that this is a current based on the perturbative calculations in terms of nuclear velocities. The transmission coefficient is decomposed into three parts: The first part T (0) (ω) is the standard transmission coefficient which instantaneously depends on the molecular geometry: with the last two terms being the dynamical corrections to the transmission coefficient due to nuclear motion. These terms have more cumbersome expressions and it would be beneficial for presentation purposes to define the following two quantities (note that we have suppressed function notation) 67 : and The coefficients A µν (x, ω) and B µ (x, ω) can then be expressed as: and where we have omitted functional dependence on x and ω in the right sides of the equations for brevity.

III. APPLICATIONS
The general theory presented in section II is now applied to several transport scenarios. The model system consists of a single molecular orbital coupled to a single nuclear degree of freedom. We structure this section as follows: first, we detail the model based expressions for viscosity, diffusion, and electric current and we also discuss the semi-analytical exact solution for the considered model. Next, we benchmark the approximate current computed within our approach against the exact results; and, finally, we consider the existence of the blowtorch effect and dynamical blockade of electric current, respectively. Note that all numerical values in the text and figures are given in atomic units.

A. Model
Before applying the theory to several model problems, we first simplify some of the general equations above in the limit of a single electronic energy level and a single nuclear DOF. We model the molecular-orbital energy as being dependent on a single classical and Cartesian DOF x(t), with the corresponding Hamiltonian being given bŷ where we have neglected electronic spin. Furthermore, we model the Hamiltonian's dependency on nuclear position in a linear fashion such that: where λ is a strength of electronic-nuclear coupling and h 0 is a fixed component of the energy level. We also choose to work in the wide-band approximation where the level-broadening functions for the left and right leads (Γ L and Γ R , respectively) are energy-independent constants. The practical result is that and hence they disappear from all expressions. Furthermore, the limit of a single electronic energy level means that all fundamental quantities -the Hamiltonian and its derivatives, the Green's functions and self-energiesbecome numbers rather than matrices and this simplifies the expressions considerably.
In these limits, the diffusion coefficient, the frictional coefficient and the adiabatic force are given by: and The coefficients A and B which determine the dynamic corrections to electric current become and This model also enables us to find the exact numerical solution of the model. As shown in the appendix, we can derive the following expression for the transmission coefficient to compute the stochastically-averaged current from the left lead according to The current equation (37) and the transmission coefficient in (36) serve as the expressions for the analytical current. Given knowledge of h(t) for a range of stochastic trajectories, then it follows that one can compute exact current trajectories for benchmarking purposes. Subsequent discussions will frequently make use of the shorthand notation J exact when referring to the analytical current.

B. Assessment of the Perturbative Expansion
We assess the performance of our theory against the exact results by benchmarking the electronic current expressions in two ways: firstly, we calculate J exact and J pert as functions of time and compare these trajectories (see figure 1); and secondly, we compute the currentvoltage characteristics using the expressions for J exact and J pert (see figure 2). In the discussions that follow, a case of good quantitative agreement between the J exact and J pert trajectories and a case of bad quantitative agreement will be considered. The reader can refer to Appendix B for a general description of the algorithm used.
In figure 1, we plot J exact and J pert as functions of time for a single Langevin trajectory once the system has reached the steady state in the case of good and bad qualitative agreement; we also plot the adiabatic current J (0) as calculated according to equations (21) and (23). Both calculations use the same set of parameters but with different masses for the classical degree of freedom: m = 5000 to show the case of good agreement, and m = 1000 for the case of bad agreement. Figure 1 a) shows that in the case of good qualitative agreement there is a perturbative current J pert that approximates closely the exact current J exact ; the exact current has been shifted slightly relative to the adiabatic current J (0) and this has been captured by the perturbative current. The case of bad qualitative agreement, on the other hand, shows that J pert poorly approximates J exact , as evidenced by the J pert trajectory drastically overshooting and undershooting the analytical J exact . For the positive voltages considered, J pert and J exact have time-averaged values that are shifted upwards slightly relative to that of the adiabatic trajectory, the J pert curve having the highest time-averaged value. The most noticeable feature of this plot is that J pert has a very large amplitude of oscillation relative to the other trajectories. Figure 2 shows a plot of the current-voltage characteristics of the device in the case of good qualitative agreement and bad qualitative agreement. Simulations are performed for the same parameters as figure 1 by considering the motion corrections to the adiabatic current J (0) as calculated according to the perturbative current J pert (as given by J pert − J (0) ) and the analytical current J exact (as given by J exact − J (0) ). The current-voltage plot is calculated in a similar fashion to the trajectories in figure 1; however, we compute these trajectories for an array of voltages and compute the mean and timeaveraged currents once the system has thermalised.
As demonstrated by figure 2, the case of good qualitative agreement exhibits much closer agreement between the time-averaged perturbative and analytical currentvoltage characteristics relative to the case of bad qualitative agreement. The time-averaged perturbative current is noticeably shifted away from the time-averaged analytical current, with this being particularly manifest in the case of bad qualitative agreement.
The characteristic frequency for nuclear motion can be computed from from spring constant and mass as ω = k/m, which ω = 0.0045 for the case of good agreement and ω = 0.01 for the case of bad agreement. Therefore, we see that, as qualitatively expected, if nuclear motion is slow compared to the electron tunnelling dynamics ω ≪ Γ, the perturbative approach produces the almost exact results, however, when the nuclear electronic timescales become comparable ω ∼ Γ the perturbative current deviates both qualitatively and quantitatively from exact values.

C. Non-Adiabatic Corrections
From our previous works we found that the presence of second-order corrections in the central-time derivatives produced a current-blockade effect in the current-voltage characteristics. Our work assumed that nuclear dynamics could be described in the equilibrium limit, meaning that we computed the average current according to the calculation where J (x, p) is the non-adiabatically corrected current computed for a given point in the phase space and is the equilibrium Boltzmann distribution function. The discussion that follows will now extend this analysis to consider classical particles that can undergo nonequilibrium motion in accordance with the Langevin equation. Our discussion will start by considering the non-equilibrium probability distributions P neq (x) that result from non-equilibrium motion and, secondly, we will analyse the resultant non-equilibrium conductancevoltage characteristics.
For this set of electronic and nuclear parameters, the conductance-voltage characteristics are calculated in the following fashion (the reader is invited to view Appendix B for a more thorough explanation): (i) for a given voltage V , n stochastic nuclear trajectories are calculated according to the Langevin equation; (ii) from these trajectories the analytical and perturbative current trajectories are calculated; (iii) once a steady-state has been reached, the remaining current trajectories are averaged; (iv) steps (i) to (iii) are repeated for values across the voltage spectrum to produce the current-voltage profile.
It will be useful in what follows to use an effective temperature quantity T ef f (x): where D(x) and ξ(x) are the previously-defined diffusion and friction coefficients. The introduction of local effective temperature quantities is not a universally valid approach and has been known to fail in many nanoscopic systems 74 . Nonetheless, we find that the effective temperature serves as a useful ratio to infer the relative strengths of the random forces and frictional forces for a given nuclear displacement x(t) and voltage V and as a predictor of steady-state characteristics. Additionally, it will be useful to use an effective adiabatic potential. This quantity describes the 'effective' classical potential that arises from the originally defined classical potential U (x) and the adiabatic force F (x) in the Langevin equation: The reader can refer to figure 3 a) and b) for a a visual representation of U ef f (x) and T ef f (x) for the parameters selected. Figure 3 c) considers two sets of probability density functions for systems with electron-vibrational couplings λ of varying strength. The first set, given by P neq (x) (the solid curves), are the non-equilibrium probability density functions; and the second set are what we will refer to as the effective probability density function P ef f (x) (the dashed curves). The density functions P neq (x) were generated by running many Langevin trajectories and building a normalised frequency distribution P neq (x) for each value of λ. The effective probability density function P ef f (x), on the other hand, is defined according to the formula: Above there is the adiabatic potential U ef f and the steady-state effective temperature T ef f relevant to each value of λ. Additionally, above we have the appropriate normalisation constant Z. The steady-state effective temperature T ef f is computed by running long trajectories and computing an average over time once a steady state has been reached; in addition, many long trajectories are run and these results are further averaged.
As can be seen in figure 3 c), the non-equilibrium distributions have non-zero expectation values x < 0 that move further to the left with increasing λ. In addition, we see that the variance of these probability distributions decreases with increasing λ. This phenomenon can be explained in large part by the adiabatic force F and its increasing prevalence with increasing λ. For the parameters set selected, the adiabatic force tilts the effective adiabatic potential to the left and, consequently, influences the classical degree of freedom to spend most of its time there. In addition, the effective temperature T ef f (x) is at a maximum at x = 0 and the particle will be influenced to move into the colder regions that are now provided by the changing effective potential U ef f (x). Another interesting point to note from figure 3 c) is the close agreement between P neq (x) and P ef f (x). This agreement has been found for a variety of parameter sets and suggests that the effective temperature T ef f (x) can be used as an effective measure of steady-state temperature characteristics.
We now turn our attention to analysing the conductance-voltage characteristics and contrasting the differences that arise in the equilibrium and nonequilibrium descriptions. As already mentioned, our previous work made use of an equilibrium description of the dynamics in accordance with a Boltzmann factor 67-69 . In this equilibrium limit, it can be shown that the transmission (43) reduces to where we emphasise that the brackets indicate averaging over nuclear velocities in accordance with the Boltzmann factor. An important consequence of this description is that terms linear in nuclear velocities and accelerations disappear once averaged and only corrections that are quadratic in nuclear velocities will impact the conductance-voltage characteristics. Figure 4 has been generated to enunciate the differences between equilibrium and non-equilibrium motion. The plots consider the equilibrium case as the blue curve and non-equilibrium motion as the orange curve; in addition, the adiabatic and non-adiabatic conductances are given by the solid and dashed lines, respectively. The first thing to notice is the difference in the adiabatic conductances between the two descriptions. These plots show in all cases that the adiabatic conductance as calculated according to the Langevin equation is less than the adiabatic conductance for the equilibrium description and, furthermore, that the discrepancy between the two increases with increasing λ. We have found that this difference is due to the position component of the nonequilibrium probability density distribution P neq (x) (as generated according to the non-equilibrium description) being broader than the Boltzmann component P eq (x) = e −U (x) kT ; additionally, we have found that the discrepancy between P neq (x) and P eq (x) increases with increasing λ. Thus, the conclusion is that the broader non-equilibrium probability density function P neq (x) in conjunction with the current surface J (0) produces a conductance profile that is lower than the equilibrium profile.
When analysing the non-adiabatic conductance of figure 4, we see that non-equilibrium nuclear dynamics have a significant effect on molecular junction electron transport properties. First, the non-equilibrium nuclear motion reduces differential conductance relative to the equilibrium case. These effects are exacerbated by changing the values of λ. Second, the second order nonadiabatic corrections act differently in equilibrium and non-equilibrium by increasing or reducing the electronic conductivity. The opposite role for non-adiabatic corrections for nonequilibrium and equilibrium nuclear dynamics can be attributed to acceleration dependent term in (22) which becomes particularly important for strongly coupled non-equilibrium electron-nuclear dynamics.

D. Landauer Blowtorch Effect and Current-Induced Bi-Stability
This section discusses an interesting phenomenon that arises in the non-equilibrium Langevin description: the natural appearance of the blowtorch effect and an important role it plays in current-induced switching between molecular configurations. The blowtorch effect is the kinetically controlled dynamics via elevation of the temperature for a portion of the potential energy surface 75,76 . As we shall now demonstrate, the Langevin description of nuclear dynamics results in steady-state nuclear temperatures that are non-constant with nuclear displacement x(t) and voltage V . These variations in steady-state temperature with displacement and voltage allow the possibility, among other things, of designing devices with voltage-dependent thermal characteristics and thermallyactivated switching. We once again find it useful for the purposes of analysis and discussion to work with the effective temperature quantity T ef f (x) defined via (40). We have provided 5 b) to see a visualisation of the effective temperature surface as a function of nuclear displacement x and voltage V . In addition, the adiabatic current J (0) has been plotted as a function of position and voltage in figure 5 c): this plot will prove to be useful when interpreting the Fano factor. This section will also make use of the effective adiabatic potential (41). Figure 5 a) shows for the parameters selected that increasing the voltage has the effect of changing the number of minima that the potential has. For a voltage of V = 0, the effective potential takes the form of a double potential as seen in figure 5 a). Increasing the voltage results in the formation of a central minima and the effective potential turns into a triple potential, an example of which is provided in figure 5 a) for a voltage of V = 0.17. Increasing the voltage further sees the effective potential morph into a single potential as the central potential dominates. Figure 5 b) shows, for the parameters selected, the relative effective temperatures as a function of position. By comparing figure 5 a) and b), we see a uniform temperature profile across the potential (the temperature is the background bath temperature in accordance with the Fluctuation-Dissipation theorem) where, upon increasing the voltage, we see that the central minima achieves relatively higher effective temperatures compared to its neighbouring wells.
We now analyse the motion of a classical degree of freedom x(t) according to a Langevin description. Figure  a) shows a nuclear trajectory that is constrained to oscillate within the previously-described system for a voltage of V = 0 and b) shows a voltage of V = 0.17. We see that the classical particle is restrained to oscillate within one of the local minima in U ef f (x) without an applied voltage: the potential barrier between the two wells and the low local temperatures does not provide enough energy to the particle to make a switch. At non-zero voltages, on the other hand, shows a particle that mimics switching behaviour between the far-left well and the far-right well. The potential barriers between wells and the local temperature provide enough energy for a particle to make a transition; note, however, how the particle mimics the switching behaviour of a double potential well as it predominantly switches between the far-left and the far-right wells. This is due to the high temperatures in the central potential relative to the low temperatures in the neighbouring wells: a particle entering the central potential will gather enough energy to make a transition much more quickly than a particle entering the far-left or the far-right wells. The corresponding current trajectory in 6 c) shows non-zero current for a voltage of V = 0.17 and 'spikes' of non-zero current. This can be explained by observing the effective adiabatic potentials in figure 5, where we see that the left and right minima correspond to regions of zero current for the parameters selected. On the other hand, nuclear positions of around x = −1 correspond to regions of non-zero current. A classical particle trapped in the left or right potential will therefore result in zero current flowing across the system and it is the switching of the classical particle among the left and right well, and its localisation to the central well, that will produce nonzero current. Notice that this 'current window' widens somewhat linearly with voltage; furthermore, notice how the magnitude of the current increases with the voltage: for low voltages the magnitude of the current is relatively  7: In a) we have the waiting times of the classical particle in the far-left well, the central well and the far-right well where waiting times are presented in terms of a ratio from 0 to 1. In b) there is the switching rates of the classical particle amongst the three wells.
sensitive to the voltage in comparison to high voltages.
In summary, switching of the classical particle amongst the wells corresponds to the switching of the electric current on and off: increasing the voltage will increase the current window (the portion of the trajectory that corresponds to non-zero current) and the magnitude of the current, and the formation of the central well will further increase the portion of the trajectory that corresponds to non-zero current.
The switching behaviour of the classical degree of freedom x(t) amongst the wells as a function of voltage is more succinctly expressed in terms of the waiting times between the three wells, as seen in figure 7. We define the waiting times as the time spent in each of the respective wells : for example, if the left and right wells have waiting times of 0.2 and 0.8 respectively, then this would indicate that, for a given trajectory, the particle spent 20% of its time in the left well and 80% of its time in the right well. Calculations of the waiting times separate the wells for a given effective potential according to the maxima that occur.
Voltages ranging from V = 0 to approximately V = 0.06 correspond to the effective adiabatic potential U ef f (x) emulating a double potential; voltages in the range of V = 0 to approximately V = 0.01 result in the nuclei remaining localised in a single well (the single well that the particle remains in corresponds to the initial condition set for the particle), while voltages ranging from V = 0.01 onwards demonstrate the particle accessing multiple wells. As we increase the voltage to a point greater than approximately V = 0.06, the effective adiabatic potential begins to look like a triple potential and the particle begins to spend an appreciable amount of time in the newly-formed central potential. The higher the voltage, the more prominent the central potential becomes, and we see the particle spends increasing amounts of time in the central potential. The horizontal line in figure 7 corresponds to a voltage of V = 0.17 and the effective adiabatic potential given by figure 5. This scenario corresponds to the triple potential having three minima of approximately equal height; note, however, that the particle spends minimal time in the central potential due to the high temperatures. We also note for increasing voltages that the time in the far-left and far-right wells is split evenly as a consequence of the symmetry of the well and its temperature profile.
We finish section III D by considering current noise. We compute the current-current autocorrelations and the Fano factor as a function of voltage. We calculate n Langevin trajectories for a given voltage. For a given trajectory i, we compute the deviation from the mean according to the formula: Note that we are using the variable J generally for any current quantity of interest: for the calculations that follow, we will take J to be the perturbative current J pert and the adiabatic current J (0) . We then compute the current noise due to nuclear motion according to the expression: where we have ignored quantum mechanical cross terms in the calculation. This in turn allows us to calculate the Fano factor F i for the i th Langevin trajectory according to: Above we have the time-averaged current J i and the Fourier transform of the mechanical noise S i (ω = 0) as evaluated for ω = 0. Depending on the requirements of the calculation, the aforementioned procedure can be replicated for all of the n Langevin trajectories to produce an averaged Fano factor: If Fano factor is smaller that 1 (F < 1)the electron transport is sub-Poissonian statistical process, if F = 1 it is Poissonian, and if F > 1 the transport is super-Poissonian.
We plot of the Fano factor as a function of voltage in figure 8 for the previously-described parameters. The plot consists of the adiabatic contribution to the Fano factor and the total Fano factor (as calculated by correcting the adiabatic Fano factor to the second-order in the central-time derivatives), as provided by the blue and orange curves, respectively. For the purposes of explanation the plot has been separated into several distinct regions, each colour coded and separated by a vertical line: the red region farthest to the left indicates the voltages for which the effective potential is a double potential; the orange region denotes voltages that generate triple potentials; and the blue region indicates when the system parameters have formed a single potential.
Considering now the adiabatic Fano factor in figure 8, we see that voltages in the range of V = 0 to approximately V = 0.01 are accompanied by Fano factors of approximately zero. The reason for this is intuitive: for the voltage range of V = 0 to approximately V = 0.01, the motion of the classical degree of freedom is localised to a single potential (see figure 7 for the waiting times) due to a shortage of energy that would be required for switching, leading to trajectories that are poorly autocorrelated. As the voltage is increased from V = 0.01 to around V = 0.04, we see that the Fano factor increases in a linear fashion. We explain this behaviour according to positive contributions of the autocorrelations outweighing the negative contributions in the following manner: as the voltage is increased from V = 0.01 to around V = 0.04, the current window is steadily increasing in width in a linear fashion with voltage and the current flowing across the junction is rapidly increasing with voltage, ultimately leading to rapidly increasing contributions; however, the switching rate is relatively low in this region (and is also seemingly increasing linearly with voltage) and negative contributions to the autocorrelations are consequently small.
Voltages of V = 0.04 to V = 0.065 display a Fano factor that remains relatively constant with voltage. This voltage region still sees the current window linearly increasing with voltage; however, the sensitivity of the current with voltage is now steadily decreasing and the switching rate is increasing more rapidly. The result is autocorrelations with increasingly dominant negative contributions that outweigh the positive contributions and flatten out the curve. Increasing the voltages further above V = 0.065 will see the double potential forming into a triple potential, the central potential forming within the current window (see a) of figure 5). The particle can now be trapped in the central potential and spends increasing amounts of time in the central potential with voltage (see figure 7). The result is the junction pro-ducing non-zero current for increasing amounts of time and the positive contributions to the autocorrelations increase. The switching rate continues to increase; although, its negative contributions are outweighed by the previously-described positive contributions and the result is an increasing Fano factor to approximately V = 0.1. Increasing the voltages beyond V = 0.1 results in switching rates that begin to dominate the autocorrelations: the negative contributions continually overpower the positive contributions and the Fano factor drops.

IV. CONCLUSIONS
In this paper we have extended our description of electronic transport in molecular junctions and its interaction with nuclear degrees of freedom. Using the NEGF formalism and a Wigner expansion to extract and separate fast and slow time-scales, we describe the electronic transport across the junction according to an extended current formula that accounts for the velocities, kinetic energies and accelerations of the classical degrees of freedom. In addition, we describe the nuclear degrees of freedom according to a set of Langevin equations with the adiabatic force, frictional coefficient and diffusion coefficients being derived self-consistently within the NEGF formalism.
As a first step in applying the theory to molecular junctions, we benchmark the extended current expressions against an analytical expression for the electric current and explore the parameter range for which the perturbative electric current expressions are in qualitative agreement. It is found that the perturbative expressions qualitatively agree with the analytical current when the ratio Ω Γ is small and disagree more severely as Ω Γ increases. The current characteristics in an equilibrium nuclear dynamics are then compared to fully non-equilibrium description. It is found that terms in the extended current expression do not vanish in the non-equilibrium description and, as a result, we observed significantly different current characteristics such as constructive contribution from inelastic scattering in the conductance characteristics relative to the destructive contribution predicted by equilibrium descriptions. In addition, we used the effective temperature extracted from electronic viscosity and diffusion coefficient to show that the probability distribution of finding the classical particle could be approximated with a Boltzmann-like probability distribution with some fictitious temperature. In the final section, we showed that the non-equilibrium description of nuclear motion can give rise to the Landauer blowtorch effect via the emergence of effective adiabatic potentials with multiple minima. The particular example that we considered consists of an effective potential that varies from a single to a triple potential with voltage and we explored the impacts that the Landauer blowtorch effect has on the trajectories, current, current switching and Fano factor.

DATA AVAILABILITY
The data that supports the findings of this study are available within the article.
Next, we employ wide-band approximation, transform self-energies to the energy domain whilst leaving molecular Green's functions time-dependent The current should satisfy the continuity equation at each time moment t where N is the total number of electrons in the molecule at time t. Then, using the continuity equation we can write left current as which yields Let us average the above equation over the Langevin trajectory. The average dN dt x = 0 (A6) obviously disappears, otherwise the molecule will accumulate or loose charge continuously. The averaged current becomes Notice, that the averaged electric current does not depend on time, since the retarded Green's function will depend on relative time only once averaged over the stochastic realisations. Our next goal is to find the explicit expression for the retarded Green's function to enter it into the expression to electric current. We begin with the equation of motion for the retarded Green's function it is reduced in the wide band approximation to This differential equation can be resolved analytically and yields Substituting retarded Green's function (A10) into the expression for current (A7) gives Changing variables of integration to τ = t − t ′ and taking into account that quantity e −i t t ′ dt1ǫ(t1) x depends on relative time only once averages over realisation of the stochastic process we arrive to

Appendix B: Numerical Algorithm Description
This section gives an outline of the technical details of numerical algorithm used in section III.
The integration of Langevin equation in the absence of fluctuation dissipation theorem and with position dependent viscosity and diffusion coefficient requires special care. We used recently proposed algorithm to integrate the Langevin equation and produce stochastic trajectory (x(t), p(t)) 77 . This process is repeated n times to produce a total of n Langevin trajectories.
The calculations of expectation values for various functions of position and momentum f (x, p) can be easily performed by the averaging over Langevin trajectories. Since the acceleration is not available directly from the used integration algorithm, the acceleration dependent term in (22) is computed as where stochastic force is omitted from the average since it will give zero contribution to the average for deltacorrelated noise.