Spectroscopy of the transition-rate matrix for molecular junctions: dynamics in the Franck-Condon regime

The quantum master equation applied to electronic transport through nanoscopic devices provides information not only on the stationary state but also on the dynamics. The dynamics is characterized by the eigenvalues of the transition-rate matrix, or generator, of the master equation. We propose to use the spectrum of these eigenvalues as a tool for the study of nanoscopic transport. We illustrate this idea by analyzing a molecular quantum dot with an electronic orbital coupled to a vibrational mode, which shows the Franck-Condon blockade if the coupling is strong. Our approach provides complementary information compared to the study of observables in the stationary state.


I. INTRODUCTION
Recent progress in nanotechnology allows to fabricate transistors with a single molecule forming the active region. [1][2][3][4] The transport properties of such devices have been studied extensively both experimentally and theoretically. Theoretical approaches have been reviewed by Andergassen et al. 5 and by Zimbovskaya and Pederson. 6 The transport properties of single-molecule devices are often dominated by strong interactions. Electronelectron interaction, which leads to the Coulomb blockade, 7,8 and the coupling of electrons to vibrational modes 7-10 have to be taken into account. We assume the nuclear motion to be slow on the time scale of electronic transitions, which is the case for most, but not for all, single-molecule devices studied so far. 11 The electron-vibron coupling is then due to the change of the equilibrium nuclear configuration with the electronic state, i.e., the Franck-Condon effect: Even if a molecule is initially in a vibrational eigenstate, after an electronic transition it will be in a superposition of vibrational eigenstates belonging to the new electronic configuration. The probability to end up in any one of them is determined by Franck-Condon matrix elements between vibrational eigenstates for the old and new electronic configurations. [12][13][14][15][16] These matrix elements can be very small if the equilibrium value of the relevant normal coordinate changes strongly with the electronic state. Consequently, the rates for electrons tunneling into or out of the molecule and thus the current can be strongly suppressed-this is the Franck-Condon blockade. [14][15][16][17][18][19][20][21] The dynamics is also unusual in this regime: Electrons tunnel in avalanches separated by quiet time intervals. 14,15 The avalanche-type transport is self-similar on intermediate timescales and also leads to a strong enhancement of the zero-frequency Fano factor. [14][15][16] Strong interactions and states far from equilibrium make the description of molecular devices difficult in general. 5,6 We here focus on the case of weak hybridization between the molecule and the leads. In this case, master-equation (ME) approaches 13,14,[16][17][18][19][22][23][24] or the equivalent real-time diagrammatic approach [25][26][27] can be employed. In principle, the ME takes into account all interactions in the molecule but requires an approximate treatment of the coupling to the leads.
The simplest non-trivial version of the ME treats the hybridization to second order (sequential tunneling) and neglects off-diagonal components of the reduced density matrix in the eigenbasis of the molecular Hamiltonian. 13,14,22 The neglect of the off-diagonal components is generally not justified if some of the eigenstates are degenerate, since then the choice of eigenstates in the degenerate subspaces is arbitrary. In the absence of spindependent terms in the Hamiltonian, spin degeneracy is always present. It is then appropriate to retain all secular components of the reduced density operator, 18,28,29 i.e., all diagonal matrix elements and off-diagonal matrix elements between degenerate states.
Going beyond sequential tunneling, at fourth order one obtains cotunneling and pair tunneling as well as a contribution to the life-time broadening of the molecular levels. 29 Koch et al. 16 include cotunneling but consider only the diagonal components of the reduced density matrix. However, integrating out the off-diagonal terms of second order generates contributions of fourth order to the transition rates between the diagonal components, which are taken into account by Leijnse and Wegewijs 18 as well as Koller et al. 29 The dynamics of charge transport through molecules, in particular in the Franck-Condon regime, have so far mostly been studied by considering the current-noise spectrum and the full counting statistics. [14][15][16]30 For the Franck-Condon regime, these studies have been augmented by real-time Monte Carlo simulations. 14-16 Recently, Donarini et al. 31 have considered a harmonically varying bias voltage. Assuming spinless electrons and employing the Markov and sequential-tunneling approximations, they find hysteretic behavior of current and charge. 31 In the present paper, we propose and illustrate a different approach to the dynamics: We study the eigenvalue spectrum of the transition-rate matrix appearing in the ME. Since these eigenvalues describe the decay rates and oscillation frequencies of deviations of the reduced density matrix from the stationary solution, they provide a complementary view of the dynamics. More specifically, we will consider the part of the spectrum that is small in absolute value, corresponding to states that decay or oscillate slowly. We use the well-studied molecular transistor with a single vibrational mode as a testbed for our idea.
The remainder of this paper is organized as follows: In Sec. II the model is briefly introduced and the ME and the relevant approximations are discussed. The transition-rate matrix is introduced and the physical interpretation of its eigenvalues is given. Results for the eigenvalue spectra in various regimes are presented and discussed in Sec. III. Section IV summarizes the main results and draws some conclusions.

A. Anderson-Holstein Hamiltonian
In the following we consider a device containing a single molecule coupled to two electrodes. The device is described by the Anderson-Holstein Hamiltonian H = H leads + H mol + H t , 14,22,32,33 where refers to the noninteracting leads. For simplicitly, we assume that both electrodes have identical band structures and a constant density of states. The operator c † νkσ creates an electron in lead ν = L, R with momentum k and spin σ. The molecule is described by where d † σ creates an electron with spin σ and energy d in the molecular orbital, is the corresponding number operator, and b † is the raising operator of a harmonic vibration mode. In a break-junction device, the on-site energy can be tuned by a gate voltage, which is absorbed into d . Finally, the tunneling between the molecule and the leads is described by where N 1 is the number of sites in one lead, which is often absorbed into t ν . 34 We here assume that the tunneling matrix elements t ν are independent of momentum and spin in the relevant energy range and that the contacts are symmetric.

B. Master equation
The reduced density operator of the molecule is where ρ(t) is the full density operator of the molecule and the leads and the trace is over all many-particle states of the leads. The full density operator satisfies the von Neumann equation There are various ways to obtain a ME starting from Eq. (5). 35 Under the condition that the system was in a product states with the leads in (separate) thermal equilibrium at some initial time t 0 , ρ(t 0 ) = ρ mol (t 0 ) ⊗ ρ 0 leads , one can derive a ME that is local in time. 27,[35][36][37][38][39][40] This so-called time-convolutionless ME reads (6) where the commutator induces the bare time evolution of the decoupled molecule and L(t, t 0 ) is a linear superoperator describing the coupling to the leads. The right-hand side of the ME is linear in ρ mol so that we can rewrite it as a product involving a superoperator A.
We write ρ mol in the basis of eigenstates |n, q to the eigenvalues E nq of H mol , where n ∈ {0, ↑, ↓, ↑↓} specifies the electronic state and q is the quantum number of the vibration. The eigenenergies are where n d (n) = 0, 1, 2 is the number of electrons in the electronic state n. The matrix elements of the reduced density operator in this basis are denoted by ρ nn qq ≡ n, q| ρ mol |n , q .
The ME (6) written in components reads L nn ,n n qq ,q q ρ n n q q .
Expressed precisely, our goal is to find the eigenvalue spectrum of the superoperator A in the ME (6). The physical interpretation of the eigenvalues becomes clear if one inserts the ansatz into Eq. (6). Here, α is a complex number and ζ α is an operator on the molecular Fock space. This leads to the eigenvalue equation The ansatz (10) indeed solves the ME if α is an eigenvalue of A. Then, Re α is the negative of the decay rate of the corresponding solution, while Im α is the angular frequency of its oscillations. Evidently, a vanishing eigenvalue α = 0 corresponds to a stationary state. The stationary density operator is thus the right eigenvector ζ 0 if we impose the normalization condition Tr ζ 0 = 1.
Since the ME (9) has to preserve the trace of ρ mol , it must satisfy for all ρ mol . Therefore, η 0 with components η nn 0,qq ≡ δ nn δ qq is a left eigenvector of A to the eigenvalue zero. This proves that at least one stationary state exists. This solution is unique if the system is ergodic in the sense that every state can be reached from every other state by a finite number of transitions. [41][42][43][44] This condition is satisfied by our model for non-zero temperature. Thus there is exactly one eigenvalue α = 0.
What is the meaning of the other right eigenvectors ζ α for α = 0? These eigenvectors are not well-formed density matrices since they have zero trace. This follows from the fact that η 0 is a left eigenvector to eigenvalue zero. Since the left and right eigenvectors to different eigenvalues are orthogonal, one has for all right eigenvectors ζ α to non-vanishing eigenvalues 0 = Tr η † 0 ζ α = nqn q η n n 0,q q * ζ nn α,qq = nq ζ nn α,qq . (13) Thus one cannot interpret ζ α as a density matrix. However, linear combinations of the form with constants c α chosen such that ρ mol is hermitian, have unit trace and statisfy the ME (6). As long as ρ mol is a positive matrix, it is a permissable time-dependent density matrix. Hence, the eigenvectors ζ α for α = 0 describe deviations of possible density matrices from the stationary state. Their time dependence is governed by the eigenvalue α. One can show that all eigenvalues α = 0 have negative real parts, 45 i.e, they describe deviations from the stationary state that decay for t → ∞.
In practice, approximations are needed to obtain the superoperator A. We employ the sequential-tunneling and secular approximations. While these are standard for the study of stationary states, we have to show that they are justified for our purpose of obtaining the eigenvalue spectrum. We assume weak coupling between molecule and leads and expand the right-hand side of Eq. (6) in powers of t L,R . For the tunnel Hamiltonian H t , only even powers of t L,R enter. One can thus write Here, is the bare time evolution from Eqs. (6) and (9). Equation (9) shows that A (0) is diagonal in the basis {|n, q n , q |}. We split the reduced density operator into a secular part ρ s and a non-secular part ρ n and expand both, This expansion allows us to write down the ME order by order in t 2 L,R . First, we discuss the stationary state, for which the left-hand side of the ME (6) vanishes. At order zero, we obtain Since A (0) is diagonal and according to Eq. (9) gives zero when acting on the secular part ρ multiplies any non-secular component of the reduced density operator by a non-vanishing factor. At second order we thus find The leading secular components are thus of order zero, while the non-secular components are at least of order two. From the ME up to second order in t L,R , one cannot obtain the second-order corrections to the secular part, ρ (2) s , this requires to go to fourth order. 18,29 The leadingorder stationary density operator is thus the solution of with the non-secular components vanishing, ρ (0) n = 0. Next, we turn to the dynamics. In the ME (9), the time derivative of all non-secular components of ρ mol picks up imaginary factors −(i/h) (E nq − E n q ). They are large compared to the typical scale introduced by the coupling since we assume the coupling to be weak and exclude accidental near-degeneracies. To put it differently, the zero-order superoperator A (0) has eigenvalues α (0),nn qq = −(i/h) (E nq − E n q ) to eigenstates |n, q n , q |. These eigenvalues are zero for secular components and large and purely imaginary for non-secular components. As long as the perturbative expansion in t L,R is justified, the full superoperator A will have eigenvalues close to the zero-order ones, thus some will be small in absolute value, while the rest will have a large imaginary part. We are interested in the part of the spectrum with small absolute value. This is the part lacking large imaginary parts from order zero. Consequently, to find the small eigenvalues, one has to consider the secular sector. To leading order, the eigenvalues are then given by the second-order superoperator A (2) . Consequently, we have to solve the ME for the secular density operator ρ s . The derivation of A (2) is standard and we omit the details. 14,16,22,27,28,46,47 Assuming that at some initial time t 0 the molecule and the leads are in a product state with the leads in separate equilibrium, ρ(t 0 ) = ρ mol (t 0 )⊗ρ 0 leads , and taking the limit t 0 → −∞, one obtains to second order in t L,R , We expand the nested commutators and take the trace over the lead degrees of freedom using (23) where ξ νk ≡ k − µ ν , µ ν is the chemical potential in lead ν, and f (ξ) is the Fermi function. The chemical potentials in the left and right leads satisfy µ R − µ L = eV . We assume that the device is symmetric, i.e., µ L = −µ R = −eV /2. In the basis of eigenstates |n, q , the ME has the form R n n q q * ρ nn qq + n n q q R nn ,n n qq ,q q ρ n n q q .
Since we are only interested in the secular sector, we can drop the first term on the right-hand side, which corresponds to A (0) . The expressions for the rates R also simplify in this sector. The rates appropriate for secular ρ mol read R nn ,n n qq ,q q where the rate Γ ≡ 2πN L,R |t L,R | 2 /h describes the coupling to the electrodes with densities of states N L,R , which are assumed to be constant, and D σ nn ≡ n|d σ |n , D †σ nn ≡ n|d † σ |n are matrix elements of the electronic operators. The Franck-Condon matrix elements F qq ≡ q|e −λ(b † −b) |q read explicitly 14,16,18,22 and In Eqs. (25) and (26), the matrix elements D σ nn and D †σ nn always appear in combinations corresponding to the creation and annihilation of electrons of the same spin σ. This leads to the vanishing of the rates for certain combinations of electronic states. Furthermore, the only off-diagonal secular components of ρ mol are ρ ↑↓ qq and ρ ↓↑ qq for all q. Thus, the rates relevant for the secular sector simplify to Note that principal-value integrals, which plague the ME for the full reduced density operator, cancel in the secular sector. Also, the rates in Eqs. (28) and (29) are real. The diagonal components of the ME (24) then simplify to while for the off-diagonal secular components we obtain The equations for the diagonal and the off-diagonal components thus decouple and the off-diagonal components exhibit simple exponential decays.
Averages of local observables such as the electron number in the molecule can be obtained directly from ρ nn qq .
The current requires a different approach. In the following, we consider a charge current flowing from left to right as positive. The operatorÎ ν for the current between the lead ν = L, R and the molecule iŝ where the numerical value of ν is +1 (−1) for the left (right) lead, −e is the charge of the electron,N ν is the number operator of electrons in lead ν, and H ν t is the part of the tunneling Hamiltonian H t involving lead ν. The current is then I ν = TrÎ ν ρ, where Tr is the trace over the full Fock space. Under the same assumptions as used for the rates above one obtains where the sum over n, q, n , q only runs over secular components.
For the numerical calculations, we cut off the ladder of harmonic-oscillator states so that 0 ≤ q ≤ q max . Then the dimension of the molecular Fock space is 4(q max + 1). The secular part of ρ mol has 6(q max + 1) componentsthe 4(q max +1) diagonal ones and 2(q max +1) off-diagonal ones of the form ρ ↑↓ qq and ρ ↓↑ qq . In the secular sector, the density operator can thus be represented by a 6(q max +1)component vector and the superoperator A (2) by a real 6(q max + 1) × 6(q max + 1) transition-rate matrix.
It is a common problem of the ME approach that the matrix A (2) can be ill-conditioned, i.e., the ratio of the largest to the smallest eigenvalue can be very large. This is here due to the Fermi functions and Franck-Condon matrix elements in Eqs. (28) and (29), which span many orders of magnitude. For this reason, black-box diagonalization routines often fail to distinguish the true stationary state from eigenvectors to very small eigenvalues. To overcome this difficulty, we use Mathematica 48 to solve the eigenvalue problem with high precision. We adapt the number of digits in the calculation such that it is larger than the L ∞ condition number of the matrix by at least 12, which should give results for the eigenvalues and eigenvectors with on the order of 12 significant digits.
Before analyzing eigenvalue spectra in the next section, we comment on their dependence on the cutoff q max . We find that adding highly excited vibrational states only adds eigenvalues to the middle and upper part of the spectrum. If q max is not too small, the spectrum at small magnitudes does not change significantly with q max . This is plausible since adding highly excited vibrational states should not introduce additional slow relaxation channels. We use q max = 30, unless noted otherwise.
It would be possible extend the analysis to higher orders in the tunneling amplitudes t L,R . This would require to solve the eigenvalue equation (11) perturbatively. Compared to time-independent perturbation theory for Hamiltonian systems, this is more complicated since the superoperator A is not hermitian. 49 In this work, we obtain the eigenvalues α up to second order. The next nonvanishing contribution is of fourth order. At this order, the restriction to the secular sector is not possible, we require A (2) for all secular and non-secular states, and we also need A (4) within the secular sector. The stationary solution at this order has been studied before. 18,29 We leave the spectral analysis for a future work.

III. RESULTS
In the present section, we present results for the eigenvalue spectrum in the regimes of the transmitting quantum dot, the Coulomb blockade, and the Franck-Condon blockade. It is shown that the spectrum differs qualitatively between these cases and is characteristic for each. We compare the information that can be gained more conventionally from observables such as current and charge in the stationary state.   Fig. 1.

A. Transmitting regime
First, we consider the transmitting molecular device. We tune the effective on-site energy d −λ 2h ω v in Eq. (7) to zero so that resonant tunneling is possible at vanishing bias voltage. Figure 1 shows the current I L , the electronic occupation number n d , and the vibron excitation q as functions of the bias voltage in the stationary state for relatively small electron-vibron coupling λ = 1. The current increases with characteristic steps. 13,14,22 There is a step at zero bias since the device is on resonance. The steps at non-zero bias result from inelastic tunneling under excitation of 1, 2, . . . vibron quanta, 13,14,22 in agreement with the observed average excitation q . The average electronic occupation changes only weakly and mostly above an energy scale on the order of U . This weak dependence is due to a small admixture of doubly occupied states at higher bias voltages.
The real and imaginary parts of the eigenvalues α of the transition-rate matrix are shown in Fig. 2 for the same parameters. The eigenvalue zero is always present, as it must be. The real and imaginary parts of the other eigenvalues reflect the step positions from Fig. 1, except for the zero-bias step. Note that at any voltage, most eigenvalues have vanishing imaginary parts. The eigenvalues with non-vanishing imaginary part form complexconjugate pairs, since the transition-rate matrix is real. The non-vanishing imaginary parts are typically small compared to the real parts. This means that the decay time of the corresponding deviations from the stationary state is much shorter than their oscillation period. Similar behavior is found for randomly distributed transition rates, where it is essentially a consequence of different scaling of the real and imaginary parts with the dimension of the molecular Fock space. 44 At zero bias, all eigenvalues are real. This has to be the case since for V = 0 the molecule is coupled to an equilibrium bath and all transition rates satisfy detailed balance. [41][42][43] This is easily confirmed by checking that the non-vanishing rates in Eq. (30) satisfy R n n,nn q q,qq /R nn ,n n qq ,q q = e (Enq−E n q )/k B T at V = 0. Then for the diagonal components, the transition-rate matrix A can be written in the form where β = 1/k B T is the inverse temperature and R 0 ij = R 0 ji . Introducing the diagonal superoperator O with components O ij = δ ij e βEi/2 , one easily sees that OAO −1 is real and symmetric and therefore has only real eigenvalues. Since this is a similarity transformation, A has the same real eigenvalues. 41,42,50 The previous argument only applies to the diagonal components. However, the off-diagonal secular components show a simple exponential decay anyway, as expressed by Eq. (31).
Concerning our goal of characterizing different regimes in terms of their eigenvalue spectra, the crucial observation is that the spectrum shows a clear gap in the real part. Thus there are no slow modes-all deviations from the stationary state decay with rates that are on the order of the characteristic rate Γ.
It is interesting to analyze the character of the stationary state and of the deviations that decay most slowly. At V = 0, the stationary (equilibrium) state is a mixture of the microstates |n, q = |0, 0 , |↑, 0 , |↓, 0 with equal probabilities, except for exponentially small thermal occupations of higher-q and doubly occupied states, see Figs. 1(b) and 1(c).
As discussed, the eigenvectors to non-vanishing eigenvalues represent deviations from the stationary state. The components ζ nn α,qq with the largest magnitudes characterize the microstates that have the largest weight in a given deviation. At V = 0, we find that the eigenvalue with the smallest non-vanishing magnitude is actually threefold degenerate-it corresponds to three linearly independent deviations. The corresponding subspace is spanned by the hermitian matrices |↑, 0 ↑, 0|−|↓, 0 ↓, 0|, |↑, 0 ↓, 0| + |↓, 0 ↑, 0| and −i |↑, 0 ↓, 0| + i |↓, 0 ↑, 0|. These matrices can be written as where σ s are the Pauli matrices and |0 0| is the projection operator onto the q = 0 vibron state. The expression (35) shows that the slowest deviations represent spin polarizations in the x, y, and z direction in the singly occupied sector. Evidently, spin polarizations decay most slowly. We will return to this point below.  At the first step in Fig. 1, where eV /2 ≈hω v , there is a crossing in the spectrum in Fig. 2(a) and we thus expect the deviation with the smallest decay rate to change in character. Figure 3 shows the secular components of the stationary density matrix and of the slowest deviation at eV = 3hω v . Compared to V = 0, the stationary state obtains finite probabilities for low-lying vibron excitions in the sectors of electronic occupation numbers 0 and 1. The deviation with the smallest decay rate is now nondegenerate and the large components ζ nn α,qq change sign when the occupation changes between 0 and 1 and also when q is increased by unity. How can we understand this? At eV = 3hω v , q can increase at most by unity in a sequential-tunneling event. Sequential tunneling also changes the occupation by ±1. The slowest deviation is dominated by an imbalance between the probabilities of the microstates |0, 0 , |↑, 1 , |↓, 1 , |0, 2 on the one hand and of |↑, 0 , |↓, 0 , |0, 1 , |↑, 2 , |↓, 2 on the other. This imbalance relaxes slowly because endothermal transitions between any microstate from one class and any microstate from the other are thermally suppressed.
At the third step at eV ≈ 6hω v , there is another crossing in Fig. 2(a). Figure 4 shows the secular components of the stationary density matrix and of the slowest deviation at eV = 11hω v . The stationary state now contains highly excited vibrons. Also, the probabilities of doubly occupied states are comparable to those of empty and singly occupied states. The slowest deviation is nondegenerate and mainly involves a transfer of weight between weakly excited vibron states with q < ∼ 5 and highly excited states with q > ∼ 5. The significance of the number of 5 becomes clear by inspecting the spectra in Fig.  2: At eV = 11hω v , q can increase by at most 5 in a sequential-tunneling event, whereas any decrease is possible. The deviation sketched in Fig. 4(b) is mainly an imbalance between vibron states that differ in q by more than 5. Such a deviation is slow to relax by endothermal sequential tunneling since the relaxation requires more than one transition. We have checked this interpretation by following the slowest deviation to higher voltages. It retains its character but the zero of ζ nn α,qq shifts to higher q (not shown). This is expected since for higher voltages larger changes of q in a single transition become possible. Figure 2(a) shows that the slowest modes become slower with increasing voltage, due to the decrease of Franck-Condon matrix elements F qq for larger |q − q |.

B. Coulomb blockade
If the molecular energy level d is detuned from resonance, there is a non-zero excitation energy between states with different occupation numbers and the current is suppressed at small bias voltages. First, we consider the case that the molecular energy level d lies below the Fermi energy of the leads at zero bias but the addition energy d + U for a second electron lies above. In this Coulomb-blockade regime, the current is suppressed by the Coulomb repulsion U . For this regime, Fig. 5 shows the current I L , the electronic occupation n d , and the vibron excitation q as functions of bias voltage in the stationary state for relatively small electron-vibron coupling λ = 1. Compared to the transmitting regime, we observe Coulomb blockade for small bias voltages (|eV | < ∼ 2hω v ), where all three observables are approximately constant and the stationary state is an equal mixture of the degenerate singly occupied ground states |↑, 0 and |↓, 0 with exponentially small corrections. When the bias voltage reaches a certain threshold, electrons can tunnel out of the molecule so that the average occupation number decreases, see Fig. 5(b), and the current sets in, see Fig.  5(a). For higher voltages, also doubly occupied states occur with significant probability and the average occupation number increases again. Similarly to the transmitting regime, the current increases in steps due to inelastic tunneling under excitation of vibrons, as seen from the increase in q in Fig. 5(c). Figure 6 shows the real and imaginary parts of the eigenvalues for the same parameters used in Fig. 5. Clearly, as the system enters the Coulomb blockade, a non-zero real eigenvalue becomes very small. This is a threefold degenerate eigenvalue corresponding to deviations of the form (35). Thus in the Coulomb blockade, the spin polarization in the singly occupied sector decays very slowly. This is easy to understand: An electron has to tunnel out of the molecule and another electron with opposite spin has to tunnel in (or vice versa) to relax the spin. But the first tunneling process is thermally suppressed by the exponentially small tail of the Fermi function. If we were to include higher orders in t L,R in the calculation, eigenvalues with exponentially small leadingorder contribution would generically obtain a contribution of order |t L,R | 4 that is not exponentially suppressed but is still small as long as the pertubative expansion in t L,R is justified. We have shown above that the same spin deviations still decay slowly, although not with exponentially suppressed rate, in the transmitting regime. Note that the inclusion of the off-diagonal secular components of ρ mol is necessary to obtain the correct spin symmetry and degeneracy of these slow modes.
Next, we turn to the two regimes where both d and d + U lie either above or below the Fermi energy of the leads. Then, the molecular orbital in the stationary state is either predominantly empty or doubly occupied, respectively. The two regimes are related to each other by a particle-hole transformation so that the transport properties are very similar. In these regimes it is the single-particle energy rather than the Coulomb interaction that suppresses sequential tunneling. We nevertheless continue to use the term "Coulomb blockade". For a predominantly empty or doubly occupied molecular orbital, the molecular spin is essentially zero and its relaxation should not be important for the dynamics, like it was in the previous case. We plot the stationary current, electronic occupation, and vibron excitation as functions of the bias voltage for d = 4.5hω v in Fig. 7. There is now a broad regime at low bias voltage where the molecular orbital is essentially empty. Singly occupied states become available above the Coulomb-blockade threshold so that a sequential-tunneling current sets in. Vibrons start to be excited at the same point since an electron tunneling out of the molecule has sufficient excess energy to excite the vibration. The corresponding eigenvalue spectra are plotted in Fig. 8. It is striking that in this case no eigenvalue becomes small right at the threshold at eV ≈ ±7hω v -the gap in the spectrum persists into the Coulomb-blockade regime. An eigenvalue approaches zero, i.e., the gap closes, only at a voltage of eV ≈ ±5hω v . For smaller voltages, even deeper in the Coulomb-blockade regime, there are additional transitions where further eigenvalues become small. Note that the stationary observables in Fig. 7 are all exponentially suppressed here.
The stationary state is of course dominated by |0, 0 throughout the blockade regime. We now analyze the deviations that become slow as the voltage is lowered. At eV ≈ ±5hω v , a non-degenerate mode becomes slow that mainly involves transfer of weight between |0, 0 and excited vibrational, and to a lesser extend electronic, states. This slow mode is sketched in Fig. 9(a). Below the voltage eV ≈ ±5hω v , the excited-state-to-excitedstate transitions from |0, q + 1 to |↑, q and |↓,q become suppressed. In particular, the only rapid decay channel of the state |0, 1 (to |↑, 0 and |↓, 0 and then to |0, 0 ) is suppressed. Therefore, the slowest deviation mostly involves transfer of weight between |0, 1 and |0, 0 .
Next, at eV ≈ ±3hω v , another non-degenerate mode becomes slow. It is sketched in Fig. 9(b). This mode involves a transfer of weight between the two lowest vibrational states |0, 0 , |0, 1 on the one hand and mainly the next state |0, 2 on the other. It becomes slow because below this voltage also the decay of |0, 2 is suppressed. Analogously, at eV ≈hω v , a further mode sketched in Fig. 9(c) becomes slow due to the suppression of the decay of |0, 3 . If we would increase the on-site energy further by means of a gate voltage, we expect more and more slow modes to appear.

C. Franck-Condon blockade
In this subsection, we turn to the signatures of Franck-Condon blockade in the spectra. Like for the transmitting regime, we tune the effective on-site energy d − λ 2h ω v to zero. Then resonant tunneling is possible at V = 0 and any suppression is due to Franck-Condon blockade. Figure 10 shows the stationary current, electronic occupation, and vibron excitation as functions of the bias voltage for intermediate electron-vibron coupling λ = 2. We choose a larger Hubbard interaction U = 12hω v since for the previously used value of U = 6hω v , the effective interaction in Eq. (7) would become attractive. The main effect of the stronger electronvibron coupling is the suppression of the zero-bias current step in Fig. 10(a).
The corresponding eigenvalue spectra are plotted in    . 11. We find a smaller gap at low bias voltage, compared to the case of λ = 1 shown in Fig. 2. At V = 0, the smallest eigenvalue is threefold degenerate and corresponds to spin imbalances of the form (35). The slowest deviations are thus the same as for the transmitting regime, but their decay rate has become even smaller. At first glance, it might be surprising that the enhancement of electron-vibron coupling leads to suppressed spin re- laxation. The reason is that in order for the spin to relax, electrons have to tunnel in and out of the molecule. At low voltage, the only available transitions are between |0, 0 on the one hand and |↑, 0 and |↓, 0 on the other. But these transitions are now suppressed by the small Franck-Condon matrix element F 00 = e −λ 2 /2 . The next eigenvalue, which is comparable in magnitude, is not degenerate and corresponds to an imbalance between the empty and singly occupied states. It becomes slow for the same reason.
Finally, we turn to the case of even stronger electronvibron coupling λ. The effective on-site energy d −λ 2h ω v is again tuned to zero. Figure 12 shows the stationary current, electronic occupation, and vibron excitation as functions of the bias voltage for strong electron-vibron coupling λ = 4 and U = 40hω v . Due to the large value of U , a larger cutoff q max = 50 is chosen here. Note the current scale in Fig. 12(a): The current is strongly reduced in magnitude for all voltages, in particular for small ones, by the Franck-Condon blockade. [14][15][16][17][18][19] In this regime, the voltage dependence of the occupation number and of the vibron excitation are also suppressed. The corresponding eigenvalue spectra are plotted in Fig. 13. The typical real and imaginary parts have become smaller and the gap is completely filled in at all bias voltages shown here. Thus there are slow modes in the whole voltage range. At least at small voltages, the character of the slowest modes is the same as for λ = 2, though: The most long-lived deviations are spin and charge imbalances, the decay of which is suppressed by small Franck-Condon matrix elements.
The inset in Fig. 12 shows details on the small real parts on a logarithmic scale. The small real parts roughly follow a log-uniform distribution for a certain range of rates. Within this range, the probability density function is approximately P (|Re α|) ∼ 1/|Re α|, i.e., it is scale-invariant. The uniform distribution of ln |Re α| is caused by the approximately exponential dependence of the Franck-Condon matrix elements F qq on q and q for q, q λ 2 . The approximate scale-invariance of the distribution of small rates implies that the dynamics of the system within a certain time window is also scale-invariant. This is consistent with the approximate self-similarity of the time-dependent transport found by Monte Carlo simulations. 14,15

IV. SUMMARY AND CONCLUSIONS
In the present paper, we have studied the eigenvalue spectrum of the transition-rate matrix in the ME for a molecular quantum dot coupled to metallic leads. The relaxational and oscillatory dynamics of deviations of the system from the stationary state are characterized by the real and imaginary parts of these eigenvalues, respectively. We have mainly considered the small eigenvalues, which describe the slow dynamics. Conceptually, this is similar to analyzing the spectrum of low-lying eigenenergies of a Hamiltonian system. We have applied this idea to a molecular transistor with an electronic orbital coupled to a vibrational mode.
The spectra differ qualitatively between a transmitting device, a molecule in the Coulomb-blockade regime, and a molecule in the Franck-Condon-blockade regime. We demonstrate that the character of deviations from the stationary state can be analyzed by considering the large components of the corresponding eigenvectors. Some of the deviations with the smallest decay rates represent non-zero spin polarizations of the molecule. They occur in groups of three degenerate modes corresponding to polarizations in the x, y, and z direction. In order to obtain these modes, all secular components of the reduced density matrix have to be included.
In the transmitting regime, the spectrum has a gap for any bias voltage, i.e., there are no slowly decaying deviations on the scale of the sequential-tunneling rate Γ. In the Coulomb-blockade regime with predominantly single occupation of the molecular orbital, the gap in the spectrum closes since relaxation of the electronic spin becomes slow. If instead the molecular orbital is predominantly empty or doubly occupied, there is no finite spin polarization and thus these slow modes do not exist. In these cases, the gap persists into the Coulomb-blockade regime. However, deep within these regimes the gap closes and more and more modes become slow at consecutive steps. These modes become slow since excited-state-to-excited-state transitions are thermally suppressed. The dynamics here contains additional information not accessible by observables in the stationary state, which show an exponentially suppressed voltage dependence.
For stronger electron-vibron couping we find that the gap becomes small even if resonant tunneling is energetically possible, since certain transition rates are suppressed by small Franck-Condon matrix elements. In the strong Franck-Condon-blockade regime, the gap closes over a broad range of bias voltages since many deviations now decay slowly. We also find an approximately scale-invariant distribution of the slowest rates, consistent with the previously observed self-similar dynamics in real time. 14,15 In the present paper, we have concentrated on the stationary and long-lived states. The spectra obtained in this paper show additional structure that we have not discussed, suggesting that much more information can be extracted from the spectra and the eigenmodes.