Charge transfer through single molecule contacts: How reliable are rate descriptions?

Background: The trend for the fabrication of electrical circuits with nanoscale dimensions has led to impressive progress in the field of molecular electronics in the last decade. However, a theoretical description of molecular contacts as the building blocks of future devices is challenging, as it has to combine the properties of Fermi liquids in the leads with charge and phonon degrees of freedom on the molecule. Outside of ab initio schemes for specific set-ups, generic models reveal the characteristics of transport processes. Particularly appealing are descriptions based on transfer rates successfully used in other contexts such as mesoscopic physics and intramolecular electron transfer. However, a detailed analysis of this scheme in comparison with numerically exact solutions is still elusive. Results: We show that a formulation in terms of transfer rates provides a quantitatively accurate description even in domains of parameter space where strictly it is expected to fail, e.g., at lower temperatures. Typically, intramolecular phonons are distributed according to a voltage driven steady state that can only roughly be captured by a thermal distribution with an effective elevated temperature (heating). An extension of a master equation for the charge–phonon complex, to effectively include the impact of off-diagonal elements of the reduced density matrix, provides very accurate solutions even for stronger electron–phonon coupling. Conclusion: Rate descriptions and master equations offer a versatile model to describe and understand charge transfer processes through molecular junctions. Such methods are computationally orders of magnitude less expensive than elaborate numerical simulations that, however, provide exact solutions as benchmarks. Adjustable parameters obtained, e.g., from ab initio calculations allow for the treatment of various realizations. Even though not as rigorously formulated as, e.g., nonequilibrium Green’s function methods, they are conceptually simpler, more flexible for extensions, and from a practical point of view provide accurate results as long as strong quantum correlations do not modify the properties of the relevant subunits substantially.


Introduction
Electrical devices on the nanoscale have received substantial interest in the last decade [1]. Impressive progress has been achieved in contacting single molecules or molecular aggregates with conducting or even superconducting metallic leads [2,3]. The objective is to exploit nonlinear transport properties of molecular junctions as the elementary units for a future molecular electronics. While the initial experiments were operated at room temperature, low temperatures down to the millikelvin range, the typical regime for devices in mesoscopic solid state physics, are also accessible (see, e.g., [4][5][6]). This allows for detailed studies of phenomena such as inelastic charge transfer due to molecular vibrations [7][8][9], voltage driven conformational changes of the molecular backbone [10], Kondo physics [11], and Andreev reflections [6], to name but a few.
These developments have been accompanied by efforts to advance theoretical approaches in order to obtain an understanding of general physical processes on the one hand and to arrive at a tool to quantitatively describe and predict experimental data. For this purpose, basically two strategies have been followed. One is based on ab initio schemes that have been successfully employed for isolated molecular structures, such as, e.g., density functional theory (DFT). Combining DFT with nonequilibrium Green's functions (NEGF) allows us to capture essential properties of junctions with specific molecular structures and geometries [2,3,12,13]. This provides insight into the electronic formations of contacted molecules and gives at least qualitatively correct results for currents and differential conductances. However, a quantitative description at the level of accuracy known from conventional mesoscopic devices still seems to be out of reach. Furthermore, these methods are not able to capture phenomena resulting from strong correlations effects, such as Kondo resonances.
Thus, an alternative route, mainly inspired by solid state methodologies, starts with simplified models that are assumed to cover the relevant physical features. The intention then is to reveal fundamental processes characteristic for molecular electronics that give a qualitative description of observations from realistic samples, but provide also the basis for a proper design of molecular junctions to exploit these processes. Information about specific molecular set-ups appears merely in the form of parameters which offer a large degree of flexibility. In general, to attack the respective many body problems, perturbative schemes have been applied, the most powerful of which are nonequilibrium Green's functions [14,15]. However, conceptually simpler, easier to implement, and often better at revealing the physics, are treatments in terms of master or rate equations. Being approximations to the NEGF frame in certain ranges of parameters space, they sometimes lack the strictness of perturbation series, but have been extensively employed for mesoscopic devices [16] and quantitatively often provide solutions of at least similar accuracy. Roughly speaking, these schemes apply as long as quantum correlations between relevant subunits of the full compound are sufficiently weak [15]. Physically, it places charge transfer through molecular contacts in the context of inelastic charge transfer through ultrasmall metallic contacts (dynamical Coulomb blockade [17]) and in the context of solvent or vibronic mediated intramolecular charge transfer (Marcus theory) [18][19][20].
While rate descriptions have been developed in a variety of formulations before [21][22][23][24][25][26][27][28], the performance of such a framework in comparison with numerically exact solutions has not yet been addressed. The reason for this is simple: A numerical method that provides numerically exact data in most ranges of parameters space (temperature, coupling strength, etc.) has only very recently been successfully implemented in the form of a diagrammatic Monte Carlo approach [29]. Path integral Monte Carlo methods have been used previously for intramolecular charge transfer in complex aggregates [18,19] in a variety of situations, including correlated [30] and externally driven transfer [31] and, of particular relevance to the present work, transfer in the presence of prominent phonon modes [32].
The goal of the present work is to study a simple yet highly nontrivial set-up, namely, a molecular contact with a single molecular level coupled to a prominent vibronic mode (phonon) which itself may or may not be embedded in a bosonic heat bath. We develop rate descriptions of various complexity, place them into the context of NEGF, and compare them with exact solutions. The essence of this study is, astonishingly enough, that rate theory provides quantitatively accurate results for mean currents over a very broad range of parameter space, even in domains where they are not expected to be reliable.

Results and Discussion
In subsection 1 we define the model and the basic ingredients for a perturbative treatment. A formulation which closely follows the P(E) theory for dynamical Coulomb blockade is discussed in subsection 2. Nonequilibrium effects in the stationary phonon distribution are analyzed in subsection 3 based on a dynamical formulation of charge and phonon degrees of freedom. The presence of a secondary bath is incorporated in subsection 4 together with an improved treatment of the molecule-lead coupling, which is exact for vanishing electron-phonon interaction. The comparison with numerically exact data and a detailed discussion is given in subsection 5.

Model
We start with the minimal model of a molecular contact consisting of a single electronic level (dot) coupled to fermionic reservoirs, where a prominent internal molecular phonon mode interacting with the excess charge is described by a harmonic degree of freedom ( Figure 1) [15,33,34]. Neglecting spin degrees of freedom the total compound Hamiltonian is thus described by (1) Here, the T k,α denote tunnel couplings between dot level and reservoir α and contains the coupling M 0 between excess charge and phonon mode. An external voltage V across the contact is applied symmetrically around the Fermi level such that ε k,α = ε 0 (k) + μ α , with the bare electronic dispersion relation ε 0 (k) and chemical potentials μ L = +eV/2, μ R = −eV/2. Below, this model will be extended further to include the embedding of the prominent mode into a large reservoir of residual molecular and/or solvent degrees of freedom acting as a heat bath. Qualitatively, since the dot occupation d † d can only take the values q = 0 or 1, the sub-unit H D + H D,Ph + H Ph describes a two state system coupled to a harmonic mode (spin-boson model [20]). Depending on the charge state of the dot the phonon mode is subject to potentials V q (x 0 ) = (m / 2)(x 0 + l 0 q) 2 . Now, the presence of the leads acts (for finite voltages) as an external driving force alternately charging (q = 1) and discharging (q = 0) the dot, thus switching alternately between V 0 and V 1 for the phonon mode. The classical energy needed to reorganize the phonon is the so-called reorganization energy . Quantum mechanically, the phonon mode may also tunnel through the energy barrier located around x 0 = −l 0 /2 separating the minima of V 0,1 .
It is convenient to work with dressed electronic states on the dot and thus to apply a polaron transformation generating the shift l 0 in the oscillator coordinate associated with a charge transfer process, i.e., (2) with momentum operator where and b 0 are creation and annihilation operators of the phonon mode, respectively. We mention in passing that complementary to the situation here, the theory of dynamical Coulomb blockade in ultrasmall metallic contacts is based on a transformation which generates a shift in momentum (charge) rather than position [17]. Now, the electron-phonon interaction is completely absorbed in the tunnel part of the Hamiltonian, thus capturing the cooperative effect of charge tunneling onto the dot and photon excitation in the molecule, i.e., with Single charge tunneling through the device can be formally and exactly captured under weak conditions (e.g., instantaneous equilibration in the leads during charge transfer) within the Meir-Wingreen formulation based on nonequilibrium Green's functions [14,15]. For the current-voltage characteristics one finds (4) with energy dependent lead self-energies ∑ α (ε) = 2π∑ k |T k,α | 2 δ(ε -ε k ) and with the Fourier transforms of the time dependent Green's functions and . Upon applying the polaron transformation (Equation 2), one has (5) where all expectation values are calculated with the full Hamiltonian (Equation 3). Of course, for T k,α → 0, the Green's functions factorize as, e.g., with the phonon correlation (6) into expectation values with respect to the dot (D) and the phonon (Ph), respectively. Any finite tunnel coupling induces correlations that in analytical treatments can only be incorporated perturbatively. There, the proper approximative scheme depends on the range of parameter space one considers. Generally speaking, there are four relevant energy scales ∑ L/R , M 0 , k B T, and ω 0 of the problem corresponding to three independent dimensionless parameters, e.g., In the following we are interested in the low temperature domain θ > 1 where thermal broadening of phonon levels is small such that discrete steps appear in the I-V characteristics. Qualitatively, seen from the dynamics of the phonon mode, two regimes can be distinguished according to the adiabaticity parameter ∑/ ω 0 = σ: For σ < 1 the phonon wave packet fulfills, on a given surface V 0 or V 1 , multiple oscillations before a charge transfer process occurs. The electron carries excess energy due to the finite voltage, and this energy may be absorbed by the phonon to promote reorganization to the new conformation (in the classical case the reorganization energy Λ). In the language of intramolecular charge transfer this scenario corresponds to the diabatic regime with well-defined surfaces V q . In the opposite regime σ > 1 charge transfer is fast such that the phonon may undergo multiple switchings between the surfaces V 0,1 . This is the adiabatic regime. In this latter range the impact of the adiabaticity on the diabatic ground state wave functions is weak for m 0 < 1 when the distance of the diabatic surfaces is small compared to the widths of the ground states. For m 0 > 1 in both regimes electron transfer is accompanied by phonon tunneling through energy barriers separating the minima of adiabatic or diabatic surfaces. The dynamics of the total compound are then determined by voltage driven, collective tunneling processes. Master equation approaches to be investigated below, rely on the assumption that both sub-units, charge degree of freedom and phonon mode preserve their bare physical properties even in the case of finite coupling m 0 . Hence, since the model (Equation 1) can be solved exactly in the limits m 0 = 0 and σ = 0 and following the above discussion, we expect them to capture the essential physics quantitatively in the domain m 0 < 1 and for all ratios σ. We note that recently the strong coupling limit including the current statistics has been addressed as well [35,36].

Rate approach I
The simplest perturbative approach considers the cooperative effect of electron tunneling and phonon excitation in terms of Fermi's golden rule for the tunneling part . For this purpose one derives transition rates for sequential transfer according to Figure 1. A straightforward calculation for energy independent self-energies ∑ L/R (wide band limit) gives the forward rate onto the dot from the left lead (8) where f β (ε) is the Fermi distribution. Inelastic tunneling associated with energy emission/absorption of phonons is captured by the Fourier transform of the phonon-phonon correlation exp[J(t)] leading to (9) with denoting the mean values for single phonon absorption (a) and emission (e). The exponentials in the prefactor contain the dimensionless reorganization energy = Λ/ ω 0 . Apparently, inelastic charge transfer includes the exchange of multiple phonon quanta according to a Poissonian distribution. Further, one has the detailed balance relation P 0 (−ε) = e −βε P 0 (ε). For vanishing phonon-electron coupling m 0 → 0 only the elastic peak survives, thus P 0 (ε) → δ(ε). We note again the close analogy to the P(E) theory for dynamical Coulomb blockade [17]. Moreover, golden rule rates for intramolecular electron transfer between donor and acceptor sites coupled to a single phonon mode are of the same form with the notable difference, of course, that in this case one has a discrete density of states for both sites [20,22]. The fundamental assumption underlying the golden rule treatment is that equilibration of the phonon mode occurs much faster than charge transfer. In the last two cases this is typically guaranteed by the presence of a macroscopic heat bath (secondary bath) strongly coupled to the prominent phonon mode. Here, the fermionic reservoirs in the leads impose phonon relaxation due to charge transfer only. Thus, for finite voltage the steady state is always a nonequilibrium state that can only roughly be described by a thermal distribution of the bare phonon system (see below). One way to remedy this problem is to introduce a phonon-secondary bath interaction as well (see below in subsection 4). The remaining transition rates easily follow due to symmetry (10) Now, summing up forward and backward events, the dot population follows from (11) with the total rate Γ tot,0 = Γ L + Γ R + Γ' L + Γ' R and the rate for transfer towards the dot Γ D = Γ L + Γ' R obtained according to Equation 8. Note that for vanishing electron-phonon coupling M 0 = 0 one has Γ tot,0 (M 0 = 0) = ∑ L + ∑ R . The steady state distribution → p dot = Γ D /Γ tot,0 is approached with relaxation rate Γ tot,0 . For a symmetric situation ∑ L = ∑ R with ε D = 0 one shows that p dot = 1/2 independent of the voltage, while asymmetric cases lead to voltage dependent stationary populations. The steady state current is given by A transparent expression is obtained for ε D = 0, namely, Despite its deficiencies mentioned above, the golden rule treatment provides already a qualitative insight into the transport characteristics. Typical results are shown in Figure 2. The I-V curves display the expected steps at . Each time the voltage eV/2 exceeds multiples of ω 0 new transport channels open associated with the excitation of one additional phonon. For higher temperatures the steps are smeared out by thermal fluctuations. The range of validity of this description follows from the fact that a factorizing assumption for the electron-phonon correlation and an instantaneous equilibration of the phonon mode after a charge transfer has been used, which means that σ < 1 and m 0 < 1. The latter constraint guarantees that conformational changes of the phonon distribution remain small.
There are now three ways to go beyond this golden rule approximation. With respect to the phonon mode, one way is to explicitly account for the nonequilibrium dynamics, another is to introduce a direct interaction with a secondary heat bath in order to induce sufficiently fast equilibration. With respect to the dot degree of freedom one can exploit the fact that for vanishing charge-phonon coupling the model can be solved exactly.

Master equation for nonequilibrated phonons
To derive an equation of motion for the combined dynamics of charge and phonon degrees of freedom, one starts from a Liouville-von Neumann equation for the full polaron transformed compound (Equation 3). Then, applying a Born-Markov type of approximation with respect to the tunnel coupling to the fermionic reservoirs, one arrives at a Redfield-type equation for the reduced density matrix of the dot-phonon system [15]. An additional rotating wave approximation (RWA) separates the dynamics of diagonal (populations) and off-diagonal (coherences) elements of the reduced density. Denoting with the probability to find q charges on the dot (here, for single charge transfer q = 0,1) and the phonon in its n-th eigenstate, one has (14) with ν 0 = 1, ν 1 = −1 and energies . The matrix elements of the phonon shift operator read (15) where 1 F 1 denotes a hypergeometric function. The underlying assumptions of this formulation require weak dot-lead coupling σ < 1 and sufficiently elevated temperatures σθ < 1 for a Markov approximation to be valid. Although we will see below when comparing low temperature results with numerically exact solutions that this seems to be only a weak constraint.
The calculation of the steady state distribution reduces to a standard matrix inversion. One can show that for a symmetric system with ε D = 0, ∑ L = ∑ R one has = . A typical example for the mean phonon number is depicted in Figure 3. The curve is well approximated by a/m 0 with a ≈ 0.7. Apparently, <n> diverges for m 0 → 0 since then and approach constants independent of the phonon number. Upon closer inspection one finds that excitation is more likely than absorption, i.e., f (n,n + 1) > f (n,n − 1), for all 0 ≤ n ≤ N 0 (m 0 ) where N 0 (m 0 ) increases with decreasing m 0 . The opposite is true for n > N 0 (m 0 ) such that in a steady state, depending on the voltage, the tendency is to have higher excited phonon states occupied by smaller couplings m 0 . In particular, for strong coupling transitions n → n + k,k ≥ 0 are blocked at small n.
A convenient strategy to include nonequilibrium effects in the phonon distribution, sometimes used in the interpretation of experimental data, is the introduction of an effective temperature T eff . This way one could return to the simpler modeling of the previous section. However, the procedure to identify is not reliable, as Figure 4 reveals. While it clearly shows the general tendency of a substantial heating of the phonon degree of freedom induced by the electron transfer, the profile of a thermal distribution strongly differs from the actual steady state distribution.  Figure 5 shows that deviations are negligible for low voltages in the regime around the first resonant step (|eV/2| < ω 0 ), where at sufficiently low temperatures only the ground state participates such that the steady state distribution coincides with the thermal one. For larger voltages deviations occur with the tendency that for smaller couplings m 0 the nonequilibrated current is always smaller than the equilibrated one (I non < I eq ), while the opposite scenario (I non > I eq ) is observed for larger m 0 . At sufficiently large voltages, one always has I non < I eq . This behavior results from the combination of two ingredients, namely, the phonon distributions and the Franck-Condon overlaps |f n,k | 2 . To see this in detail, let us consider a fixed voltage. Then, on the one hand, for smaller m 0 the steady state distribution is broad (cf. Figure 3), such that, due to normalization, less weight is carried by lower lying states compared to a thermal distribution at low temperatures; on the other hand, for m 0 < 1 the overlaps |f n,k | 2 favor contributions from low lying states in the current (16), which is thus smaller than I eq . For increasing electron-phonon coupling m 0 > 1, the overlaps |f n,k | 2 tend to include broader ranges of phonon states also covered by , compared to those of low temperature thermal states. A voltage dependence arises since with increasing voltage higher lying phonon states participate in the dynamics supporting the scenario for smaller couplings. Interestingly, as already noted in [23] the overlaps |f n,k | 2 may vanish for certain combinations of n,m depending on m 0 due to interferences of phonon eigenfunctions localized on different diabatic surfaces V q , where q = 0,1.

Rate approach II
The assumption of a thermally distributed phonon degree of freedom during the transport can be physically justified only if this mode interacts directly and sufficiently strongly with an additional heat bath (secondary bath) realized, e.g., by residual molecular modes. Here we will generalize the formulation of subsection 2 to a situation where the secondary bath is characterized by Gaussian fluctuations. Its corresponding modes can thus effectively be represented by a quasi-continuum of harmonic oscillators for which the phonon correlation function (Equation 6) can be calculated easily (17) Here the spectral density I(ω) now describes the combined distribution of the prominent mode and its secondary bath. It is thus proportional to the imaginary part of the dynamic susceptibility of a damped harmonic oscillator [20]. For a purely ohmic distribution of bath modes, one has (18) where γ denotes the coupling between phonon mode and bath. The Fourier transform of exp(J) reads at finite temperatures (19) with the frequency Ω given by Ω = ω 0 ξ + iγ/2 and where the parameter ξ is . Further, ρ γ,e (Ω) = −ρ γ,e (Ω * ) (* means complex conjugation) and ρ γ = [ρ γ,a + ρ γ,e ]/2. In the above expression, contributions from the Matsubara frequencies in Equation 17 have been neglected, since they are only relevant in the regime γ β >> 2π, which is not studied here. Apparently, the coupling to the bosonic bath effectively induces a broadening of the dot levels γ(k + l)/2 compared to the purely elastic case (Equation 9). In the low temperature regime, where for equilibrated phonons absorption (related to k) is negligible, the widths grow proportionally to l. The presence of the secondary bath drives the prominent phonon mode towards thermal equilibrium with a rate proportional to this broadening. Hence, if the time scale for thermal relaxation is sufficiently smaller than the time scale for charge transfer, i.e., 1/τ l ≡ (∑ L + ∑ R )/γ << 1, the assumption of an equilibrated phonon mode is justified and the golden rule formulation (Equation 13) can be used with P 0 (ε) → P γ (ε). However, this argument no longer applies in the overdamped situation γ/ω 0 >> 1, where the phonon mode exhibits a sluggish thermalization on the time scale , which may easily exceed τ l .
As already mentioned above, for vanishing charge-phonon coupling m 0 = 0, the model (Equation 3) can be solved exactly for all orders in the lead-dot coupling [15]. In the frame of a rate description, one observes that in this limit the dot population (Equation 11) decays proportionally to (∑ L + ∑ R ). The golden rule version of the theory neglects this broadening in Equation 13 since it is associated with higher order contributions to the current (Equation 13). Now, recalling that P 0 (ε) reduces to a delta function for m 0 → 0, this finite lifetime of the electronic dot level is included for all orders by performing the time integral in the Fourier transform with ε → ε − i(∑ L + ∑ R )/ 2 ≡ ε − iΓ tot (M 0 = 0)/2 [see Equation 11]. In fact, this way one reproduces the exact solution (one electronic level coupled to leads with energy independent couplings), i.e., its exact spectral function. To be specific, let us restrict ourselves for the remainder of this discussion to the symmetric situation ∑ L = ∑ R ≡ ∑/2, and ε D = 0. Then, in the presence of the phonon mode (m 0 ≠ 0) the corresponding function follows from Equation 9 by replacing the delta function by i/[ε + ω(k − l) + i∑/2]. Again following the idea of a rate treatment, an improved version of this result accounting for higher order electron-phonon correlations is obtained through the decay rate Γ tot,0 (M 0 ≠ 0), instead of the bare dot level width ∑/ ≡ Γ tot,0 (M 0 = 0). Equivalently, one replaces i/[ε + ω(k − l) + i∑/ 2] → i/[ε + ω(k − l) + iΓ tot,0 /2] to arrive at an improved . We note that within a Green's function approach, and upon approximating the corresponding equations of motion, a similar result has been found in [15,33], with the difference though that instead of Γ tot,0 an imaginary part of a phonon statedependent self-energy appears. One can show that the Γ tot,0 appearing here within a rate scheme is related to a thermally averaged .
Now, an additional secondary bath can be introduced as above by combining Equation 19 with , leading eventually to (20) The width of the electronic dot level is thus voltage dependent and approaches the bare width from below for large voltages, that is lim V→∞ Γ tot,0 (V) = ∑/ . The range of validity of this scheme is the following: It applies to all couplings σ in the domain where the electron-phonon coupling is weak m 0 < 1. In particular, second order processes in σ capture cotunneling processes. For m 0 > 1 charge transfer is strongly suppressed and the phonon dynamics still occurs on diabatic surfaces for σ << 1 so that we expect the approach to cover this range as well.

Comparison with numerically exact results
A numerically exact treatment of the nonequilibrium dynamics of the model considered here is a formidable task. The number of formulations which allow simulations in nonperturbative ranges of parameter space is very limited. Among them is a recently developed diagrammatic Monte Carlo approach (diagMC) based on a numerical evaluation of the full Dyson series, which, in contrast to numerical renormalization group (NRG) methods [37], covers the full temperature range. For calculations of single charge transfer, results have been obtained with and without the presence of a secondary bath interacting with the dot phonon mode.
We note that computationally these simulations are very demanding as for each parameter set and a given voltage the stationary current for the I-V curve needs to be extracted from the saturated value of the time dependent current I(t) for longer times. Typical simulation times are on the order of several days to weeks, depending on the parameter range. In contrast, rate treatments require minimal computational effort and can be done within minutes. Here, we compare numerically exact findings with those gained from the various types of rate/master equations discussed above.
We start with the scenario where the coupling to a secondary bath is dropped (γ = 0) to reveal the impact of nonequilibrium effects in the phonon mode. The formulation for an equilibrated phonon is based on Equation 13 with P 0 replaced by in Equation 20, while the steady state phonon distribution is obtained from the stationary solutions to Equation 14. In the latter approach the intrinsic broadening of the dot electronic level due to coupling to the lead is introduced in the following way: One first determines via Equation 14 a steady state distribution . This result is used for an effective self-energy contribution (total decay rate) for nonequilibrated phonons, i.e., (21) where . We note in passing that lim V→∞ Γ tot,neq (V) = ∑/ ≡ Γ tot (M 0 = 0). Subsequently, an im-  proved result for the steady state phonon distribution at a given voltage is evaluated working again with Equation 14, but using the replacement: (22) Of course, for ∑ → 0 the standard Fermi distribution is regained. The corresponding steady state phonon distribution eventually provides the current according to Equation 16 with the same replacement (Equation 22) in this expression. The procedure relies on weak electron-phonon coupling m 0 < 1 and in principle also requires sufficiently elevated temperatures. Results are shown in Figure 6 together with corresponding diagMC data for various coupling strengths m 0 . Interestingly, the equilibrated model describes the exact data very accurately from weak up to moderate electron-phonon coupling m 0 ≈ 1, while deviations appear for stronger couplings and voltages beyond the first plateau eV > 2 ω 0 . For m 0 > 1 nonequilibrium effects are stronger and the corresponding master equation (Equation 14) gives a better description of higher order resonant steps. Moreover, as already addressed above, even in this low temperature domain the approximate description provides quantitatively reliable results.
In Figure 7 the frequency of the phonon mode is fixed and only the electron-phonon coupling is tuned over a wider range. For strong coupling (here m 0 = 2) the equilibrated (nonequilibrated) model predicts a smaller (larger) current than the exact one in contrast to the situation for smaller m 0 . This phenomenon directly results from what has been said above in subsection 3: For stronger coupling the Franck-Condon overlaps favor higher lying phonon states that are suppressed by a thermal distribution.
After all, the approximate models give not only a qualitatively correct picture of the exact I-V curves, but even provide a reasonable quantitative description in this low temperature domain.
In a next step the coupling to a secondary bath is turned on (γ ≠ 0) enforcing equilibration of the phonon mode, see Equation 20. The expectation is that in this case departures from the equilibrated model are reduced. In Figure 8 data are shown for a ratio m 0 = 4/5, where deviations occur at larger voltages, as observed in the previous figures. Obviously, due to the damping of the phonon mode the resonant steps are smeared out with increasing γ. However, the approximate model predicts this  effect to be more pronounced as compared to the exact data, particularly for stronger coupling γ/∑ > 1, while still γ/ω 0 < 1. In fact, in the limit of very large coupling only the k = l = 0 contribution to Equation 20 survives, such that at zero temperature one arrives at (23) with the current at large voltages I ∞ = e∑/4 and Γ tot,0 (V) ≤ ∑/ , where equality is approached for V → ∞. It seems that a broadened equilibrium distribution of the phonon, induced by the secondary bath according to Equation 20, overestimates the broadening of individual levels. Since the approach is exact in the limit m 0 → 0, the deviations appearing in Figure 8 are due to intimate electron-phonon/secondary bath correlations not captured by the rate approach. In the overdamped regime, i.e., γ/ω 0 > 1, the dynamics of the phonon mode slow down and may become almost static on the time scale of the charge transfer. In this adiabatic regime an extended version of the master equation (Equation 14) is not trivial since the conventional eigenstate representation becomes meaningless. It would be better then to switch to phase-space coordinates and develop a formulation based on a Fokker-Planck or Smoluchowski equation for the phonon. This will be the subject of future research.
The essence of this comparison is that, as anticipated from physical arguments already in subsection 1, a rate description does indeed provide quantitatively accurate results in the regime of weak to moderate electron-phonon coupling m 0 < 1 and for all σ. Deviations that occur for larger values of m 0 can partially be explained by nonequilibrium distributions in the phonon distribution, where, however, the master equation approach seems to overestimate this effect. In order to obtain some insight into the nature of this deficiency, a minimal approach consists of extending Equation 14 with Equation 22, a mechanism that enforces relaxation to thermal equilibrium with a single rate constant Γ 0 that serves as a fitting parameter. Accordingly, the respective time evolution equation for receives an additional term −Γ 0 [ − ] with the Boltzmann distribution for the bare phonon degree of freedom . Corresponding results for the same parameter range as in Figure 6 and Figure 7 are shown in Figure 9 and Figure 10 including comparison with the exact diagMC data. There, the same equilibration rate Γ 0 /∑ = 0.25 is used for all parameter sets. Astonishingly, this procedure provides excellent agreement over the full voltage range. It improves results particularly in the range of moderate to stronger electron phonon coupling, but has only minor impact for m 0 < 1. The indication is thus that electron-phonon correlations neglected in the original form of the master equation have effectively the tendency to support faster thermalization of the phonon. Indeed, preliminary results with a generalized master equation, where the coupling between diagonal (populations) and off-diagonal (coherences) elements of the reduced charge-phonon density matrix is retained (no RWA approximation), indicate that this coupling leads to an enhanced phonon-lead interaction and thus to enhanced phonon equilibration.