Vibration-mediated resonant tunneling and shot noise through a molecular quantum dot

# Vibration-mediated resonant tunneling and shot noise through a molecular quantum dot

## Abstract

Motivated by a recent experiment on nonlinear tunneling in a suspended Carbon nanotube connected to two normal electrodes [S. Sapmaz, et al., Phys. Rev. Lett. 96, 26801 (2006)], we investigate nonequilibrium vibration-mediated sequential tunneling through a molecular quantum dot with two electronic orbitals asymmetrically coupled to two electrodes and strongly interacting with an internal vibrational mode, which is itself weakly coupled to a dissipative phonon bath. For this purpose, we establish rate equations using a generic quantum Langevin equation approach. Based on these equations, we study in detail the current-voltage characteristics and zero-frequency shot noise, paying special attention to the advanced or postponed of the appearance of negative differential conductance and super-Poissonian current noise resulting from electron-phonon-coupling induced selective unidirectional cascades of single-electron transitions.

###### pacs:
85.65.+h, 71.38.-k, 73.23.Hk, 73.63.Kv, 03.65.Yz
1

## I Introduction

Recent progress in nanotechnology has facilitated the fabrication of single-electron tunneling devices using organic molecules. In turn, this has given rise to a large body of experimental(2); (1); (3); (4); (5); (6); (7) and theoretical(8); (9); (10); (11); (12); (13); (14); (15); (16); (17); (18) work concerning phonon-mediated resonant tunneling through a quantum dot (QD) with strong coupling to an internal vibrational (phonon) mode (IVM). In particular, Carbon nanotubes (CNT) have recently become the focus of much research interest because electronic transport measurements show that they demonstrate perfect signatures of phonon-mediated tunneling, e.g. stepwise structures in the current-voltage characteristics having equal widths in voltage and gradual height reduction by the Franck-Condon (FC) factor.(6); (7)

More interestingly, the experimental measurement of S. Sapmaz, et al.(7) has revealed some more subtle transport features in a suspended CNT connected to two electrodes, particularly the appearance of a weak negative differential conductance (NDC) at the onset of each phonon step followed by a sudden suppression of current (a strong NDC) at a finite bias voltage after several steps for a long CNT sample. Theoretically, the reason for the weak NDC is quite clear in that it is ascribed to the combined effect of strong electron-phonon coupling (EPC) and low relaxation, i.e. an unequilibrated phonon (hot phonon), and strongly asymmetric tunnel-couplings to the left and right electrodes.(8); (14); (15) Furthermore, McCarthy, et al. have theoretically predicted the catastrophic current decrease, but their calculations are based on the presumption that the EPC is dependent on the applied bias voltage. It is appropriate to explore yet another explanation of the origin of the strong NDC without such an assumption.

Up to now, most theoretical works have focused on the case of a single level coupling to the phonon mode. However, it is believed that this catastrophic current decrease is intimately related to the tunneling processes involving two electronic energy levels. Nowack and Wegewijs have considered vibration-mediated tunneling through a two-level QD with asymmetric couplings to an IVM.(16) Albeit their results display strong NDC through a competition between different FC tailored tunneling processes of the two levels, their calculations do not fully resolve the overall features of the experimental data in Ref. (7). Following a suggestion by Hettler, et al.,(19) that the interplay of strong Coulomb blockade and asymmetric tunnel-couplings of two levels to electrodes can lead to a strong NDC under certain conditions without EPC, we employ a model with two electronic orbitals having asymmetric couplings to the electrodes and with both strongly interacting with an IVM. We establish approximate rate equations at high temperature to describe resonant tunneling incorporating the unequilibrated phonon effect. Our results are in good qualitative agreement with the experiment data of Ref. (7): in the weak bias voltage region where the first molecular orbital (MO) is dominant in electronic tunneling, the current shows the FC-type steplike structure and weak NDC at the onset of each step; while with increasing bias voltage the second MO becomes dominant in tunneling, inducing a sudden rapid reduction of current whether there is EPC or not. Furthermore, our studies predict that the onset of strong NDC as a function of voltage may differ from the value at which the second MO becomes active in resonant tunneling as described in Ref. (19) due to the unequilibrated phonon effect. We ascribe this to the EPC induced selective cascades of single-electron transitions between different charge states. In addition, we also analyze the current noise properties of this system. We find that the occurrence of weak and strong NDCs is accompanied by corresponding weak and strong enhancements of zero-frequency shot noise, respectively.

To emphasize the external-voltage-driven unequilibrated phonon effect in electronic tunneling in a CNT, we consider the IVM in our model coupled to an external dissipative environment (a phonon bath),(18) and we incorporate the dissipation mechanism of the unequilibrated phonon into the ensuing rate equations on a microscopic basis. For extremely strong dissipation, our results naturally reduce to those of equilibrated phonon-mediated tunneling. Therefore, our rate equations provide a valuable theoretical framework for analysis of the mechanism by which finite phonon relaxation influences the electronic transport properties of CNT.

The outline of the paper is as follows. In Sec. II, we describe the model system that we study, and rewrite the model Hamiltonian in terms of an electron-phonon direct product (EPDP) state representation, which is suitable for the ensuing theoretical derivation. In Sec. III, we derive a set of rate equations using a generic quantum Langevin equation approach with a Markovian approximation, which facilitates investigation of uneqilibrated phonon and phonon dissipation effects on phonon-mediated electronic tunneling. In this section, we also derive the current formula using linear-response theory. In Sec. IV, we describe MacDonald’s formula for calculating zero-frequency shot noise by rewriting the ensuing rate equations in a number-resolved form. Then, we investigate in detail the vibration-mediated transport and shot noise properties of a molecular QD with two MOs in Sec. V. Finally, a brief summary is given in Sec. VI.

## Ii Model Hamiltonian

In this paper, we consider a generic model for a molecular QD (in particular the suspended CNT) with two spinless levels, one as the highest-occupied MO (HOMO) and the other as the lowest-unoccupied MO (LUMO) , coupled to two electrodes left (L) and right (R), and also linearly coupled to an IVM of the molecule having frequency with respective coupling strengths and . We suppose that this single phonon mode is coupled to a dissipative environment represented as a set of independent harmonic oscillators (phonon-bath). The model Hamiltonian is

 H=Hleads+Hmol+HB+HI, (1a) with Hleads = ∑η,kεηkc†ηkcηk, (1b) Hmol = ∑j=1,2εjc†jcj+Un1n2 (1d) +ω0a†a+∑j=1,2λjc†jcj(a†+a), HB = ∑pωpb†pbp, (1e) HI = HT+HvB, (1f) HT = ∑η,k,j(Vηjc†ηkcj+H.c.), (1g) HvB = (a†+a)∑pκp(b†p+bp), (1h)

where () is the creation (annihilation) operator of an electron with momentum , and energy in lead (), and () is the corresponding operator for a spinless electron in the th level of the QD (). denotes interdot Coulomb interaction and is the electron number operator in level . () and () are phonon creation (annihilation) operators for the IVM and phonon-bath (energy quanta , ), respectively. represents the coupling constant between electron in dot and the IVM; is the coupling strength between the IVM and phonon bath; describes the tunnel-coupling between electron level and lead . Here, we denote the density of states of the phonon bath with respect to frequency by , and define the corresponding spectral density as

 JB(ωp)=κ2pD(ωp). (2)

In the literature, the following form of the spectral density is usually considered:(20)

 JB(ωp)=κ0ωαpθ(ωp), (3)

in which is the Heaviside step function. The bath is said to be Ohmic, sub-Ohmic, and super-Ohmic if , , and , respectively. We use units with throughout the paper. In addition, it is worth noting that we do not consider the EPC-induced coupling between two MO’s in the Hamiltonian, Eq. (1d), because it is much weaker than the coupling of a single MO.

It is well-known that, in the strong electron-phonon interaction problem, it is very convenient to introduce a standard canonical transformation(21) to the Hamiltonian Eq. (1a), , with (), which leads to a transformed Hamiltonian

 ˜Hmol = ∑j˜εjc†jcj+˜Un1n2+ω0a†a, (4a) ˜HI = ˜HT+˜HvB, (4b) ˜HT = ∑η,k,j(Vηjc†ηkcjXj+H.c.), (4c) ˜HvB = (a†+a−2∑jgjnj)∑pκp(b†p+bp), (4d) with ˜εj=εj−λ2jω0 and ˜U=U−2λ1λ2ω0. Importantly, the X-operator describes the phonon renormalization of dot-lead tunneling, Xj=egj(a−a†). (4e)

Obviously, the model described by the above Hamiltonian Eq. (1a) involves a many-body problem with phonon generation and annihilation when an electron tunnels through the central region. Therefore, one can expand the electron states in the dot in terms of direct product states composed of single-electron states and -phonon Fock states. In this two-level QD system, there are a total of four possible electronic states: (1) the two levels are both empty, , and its energy is zero; (2) the HOMO is singly occupied by an electron, , and its energy is ; (3) the LUMO is singly occupied, , and its energy is ; and (4) both orbitals are occupied, , and its energy is . Of course, if the interdot Coulomb repulsion is assumed to be infinite, the double-occupation is prohibited. For the sake of convenience, we assign these Dirac bracket structures as dyadic kets (bras):(22) namely, the slave-boson kets , , and pseudo-fermion kets , . Correspondingly, the electron operator and phonon operator can be written in terms of such ket-bra dyadics as:

 cj = ∞∑n=0(e†nfjn+s¯jf†¯jndn), (5) a = ∞∑n=0√n+1(e†nen+1+∑jf†jnfjn+1+d†ndn+1), (6)

with and ( is due to anti-commutation relation of the Fermion operator). With these direct product states and dyadics considered as the basis, the density-matrix elements may be expressed as , , and .

The transformed Hamiltonian can be replaced by the following form in the auxiliary particle representation:

 ˜Hmol = ∑n[nω0e†nen+∑j(˜εj+nω0)f†jnfjn (7b) +(˜ε1+˜ε2+˜U+nω0)d†ndn], ˜HT = ∑η,k,j,n[Vηjc†ηk(e†nfjn+s¯jf†¯jndn)Xj+H.c.], (7d) ˜HvB = ∑pκp(b†p+b−p)[∑n√n+1(e†nen+1 (7g) +∑jf†jnfjn+1+d†ndn+1+H.c.) −2∑j,ngj(f†jnfjn+d†ndn)].

Based on this transformed Hamiltonian, we derive a rate equation for description of the dynamics of the reduced density matrix of the combined electron and IVM system in the sequential tunneling regime.

## Iii Quantum Langevin equation approach and rate equations

In this analysis, we employ a generic quantum Langevin equation approach,(23); (24); (25); (26); (27); (28); (29) starting from the Heisenberg equations of motion (EOMs) for the density-matrix operators , (), and :

 i˙^ρn00 = [e†nen,˜H]−=∑η,k,j,m[Vηjc†ηke†nfjmXj,nm (8c) −f†jmenX†j,nmcηk]+∑pκp(b†p+b−p) ×(√n+1e†nen+1+√ne†nen−1−H.c.), i˙^ρn11 = [f†1nf1n,˜H]−=∑η,k,m[c†ηk(Vη2f†1ndmX2,nm (8h) −Vη1e†mf1nX1,mn)+(Vη1f†1nemX†1,mn −Vη2d†mf1nX†2,nm)cηk]+∑pκp(b†p+b−p) ×(√n+1f†1nf1n+1+√nf†1nf1n−1−H.c.), i˙^ρn22 = [f†2nf2n,˜H]−=∑η,k,m[−c†ηk(Vη1f†2ndmX1,nm (8m) +Vη2e†mf2nX2,mn)+(Vη1f†2nemX†2,mn +Vη2d†mf2nX†1,nm)cηk]+∑pκp(b†p+b−p) ×(√n+1f†2nf2n+1+√nf†2nf2n−1−H.c.), i˙^ρndd = [d†ndn,˜H]−=∑η,k,m[c†ηk(Vη1f†2mdnX1,mn (8q) −Vη2f†1mdnX2,mn)−(Vη1d†nf2mX†1,mn −Vη2d†nf1mX†2,mn)cηk]+∑pκp(b†p+b−p) ×(√n+1d†ndn+1+√nd†ndn−1−H.c.), where Xj,nm = ⟨n|Xj(t)|m⟩, (8r) X†j,nm = ⟨m|X†j(t)|n⟩. (8s)

These matrix elements can be calculated as:(21)

 Xj,nm = X†j,nm (9) = ⎧⎪ ⎪⎨⎪ ⎪⎩e−g2j/2gm−nj√n!m!Lm−nn(g2j),n≤m,e−g2j/2(−gj)n−m√m!n!Ln−mm(g2j),n>m, (10)

where is the generalized Laguerre polynomial. The rate equations are obtained by taking statistical expectation values of the EOMs, Eqs. (8c)-(8q), which clearly involve the statistical averaging of products of one reservoir (phonon-bath) variable and one device variable, such as and .

To determine the products, and , we proceed by deriving EOMs for the system, phonon-bath, and reservoir operators, , , , and :

 i˙F1,nm = [e†nf1m,˜H]−=[˜ε1+(m−n)ω0]F1,nm (11b) +[F1,nm,˜HT]−+[F1,nm,˜HvB]−, i˙^ρn,n+100 = [e†nen+1,˜H]−=ω0^ρn,n+100+[^ρn,n+100,˜HT]− (11d) +[^ρn,n+100,˜HvB]−, i˙bp = [bp,˜H]−=ωpbp+[bp,˜HvB]−, (11e) i˙cηk = [cηk,˜H]−=ϵηkcηk+[cηk,˜HT]−. (11f)

The EOM for is easily obtained by Hermitian conjugation of the equations for . Formally integrating these equations, (11b)-(11f), from initial time to we obtain

 F1,nm(t) = e−i[˜ε1+(m−n)ω0]tF1,nm(0) (12d) −i∫t0dt′e−i[˜ε1+(m−n)ω0]τ[F1,nm(t′),˜HT(t′)]− −i∫t0dt′e−i[˜ε1+(m−n)ω0]τ[F1,nm(t′),˜HvB(t′)]−, ^ρn,n+100(t) = e−iω0t^ρn,n+100(0) (12h) −i∫t0dt′e−iω0τ[^ρn,n+100(t′),˜HT(t′)]− −i∫t0dt′e−iω0τ[^ρn,n+100(t′),˜HvB(t′)]−, bp(t) = e−iωptbp(0)−i∫t0dt′e−iωpτ (12j) ×[bp(t′),˜HvB(t′)]−, cηk(t) = e−iϵηktcηk(0)−i∫t0dt′e−iϵηkτ (12l) ×[cηk(t′),˜HT(t′)]−,

with . In the absence of tunnel-coupling, , we have

 Fo1,nm(t) = e−i[˜ε1+(m−n)ω0]τFo1,nm(t′), (13a) ^ρn,n+100,o(t) = e−iω0τ^ρn,n+100,o(t′), (13b) bop(t) = e−iωpτbop(t′), (13c) coηk(t) = e−iϵηkτcoηk(t′). (13d)

A standard assumption in the derivation of a quantum Langevin equation is that the time scale of decay processes is much slower than that of free evolution, which is reasonable in the weak-tunneling approximation. This bespeaks a dichotomy of time-developments of the involved operators into a rapidly-varying (free) part and a slowly-varying (dissipative) part. Focusing attention on the slowly-varying decay processes, and noting that the infinitude of macroscopic bath variables barely senses reaction from weak interaction with the QD, it is appropriate to substitute the time-dependent decoupled reservoir, phonon-bath, and QD operators of Eqs. (13a)-(13d) into the integrals on the right of Eqs. (12d)-(12l). This yields approximate results for the reservoir operators as:(23); (24); (28); (29)

 cηk(t)=coηk(t)+crTηk(t), (14a) with crTηk(t)=−i∫t0dτ[coηk(t),˜HoT(t′)]−, (14b)

where is composed of the operators in which are replaced by their decoupled counterparts (interaction picture). In fact, this is just the operator formulation of linear response theory. Similarly, the approximate results for the QD and phonon-bath are also divided into two parts:

 F1,nm(t) = Fo1,nm(t)+FrT1,nm(t)+FrvB1,nm(t) (15a) = Fo1,nm(t)−i∫t0dτ[Fo1,nm(t),˜HoT(t′)]− (15c) −i∫t0dτ[Fo1,nm(t),˜HovB(t′)]−, ^ρn,n+100(t) = ^ρn,n+100,o(t)+^ρn,n+100,rT(t)+^ρn,n+100,rvB(t) (15d) = ^ρn,n+100,o(t)−i∫t0dτ[^ρn,n+100,o(t),˜HoT(t′)]− (15f) −i∫t0dτ[^ρn,n+100,o(t),˜HovB(t′)]−, bp(t) = bop(t)+brvBp(t) (15g) = bop(t)−i∫t0dτ[bop(t),˜HovB(t′)]−. (15h)

Note that we use the super(sub)scripts and denote the reactions from tunnel-coupling and environmental dissipation, respectively.

Employing the approximate solutions of Eqs. (14a), (14b) and (15c), we can evaluate as

 I = ∑η,kVη1⟨[co†ηk(t)+cr†ηk(t)] (17) ×[Fo1,nm(t)+Fr1,nm(t)]X1,nm⟩ ≃ ∑η,kVη1⟨[co†ηk(t)Fr1,nm(t)+cr†ηk(t)Fo1,nm(t)]X1,nm⟩. (18)

The statistical averages involved here can be taken separately in regard to the electron ensembles of the reservoirs (many degrees of freedom) and in regard to the few degrees of freedom of the QD EPDP states. Accordingly, the statistical average of the product of one reservoir operator and one system operator factorizes in the averaging procedure. Therefore, the statistical average of vanishes. Moreover, in the sequential picture of resonant tunneling, the tunneling rates and current are proportional to second-order tunnel-coupling matrix elements. We thus neglect the term as it is proportional to the third-order tunnel-coupling matrix element, . The other interaction terms arise from tunneling reaction, and are of second-order of . After some lengthy but straightforward algebraic calculations, we obtain

 I = −i∑ηΓη1∫dϵ∫t0dτei(ϵ−˜ε1)τ{fη(ϵ)ρn00 (20b) −[1−fη(ϵ)]ρm11}¯X†1,nmX1,nm, in which Γηj=2πϱη|Vηj|2 (20c) denotes the tunneling strength between the molecular orbital j and lead η (ϱη is the density of states of lead η), fη(ϵ) is the Fermi-distribution function of lead η with temperature T, and ¯Xj,nm = ⟨n|Xj(t′)|m⟩, (20d) ¯X†j,nm = ⟨m|X†j(t′)|n⟩. (20e) Considering a(t)=e−iω0τa(t′), we have ¯Xj,nm = ei(m−n)ω0τXj,nm, (20f) ¯X†j,nm = e−i(m−n)ω0τXj,nm. (20g)

In the derivation of Eq. (20b), we assume that states with different phonon-numbers are completely decoherent, and , owing to big energy difference between the two direct product states and if . Moreover, a Markovian approximation will be adopted by making the replacement

 ∫t−∞dτ⟹∫∞−∞dτ, (21)

in the statistical averaging, since we are interested in the long time scale behavior of these density matrix elements.

Furthermore, can be evaluated using Eqs. (15f), (15h) with as

 J ≃ ∑pκp√n+1⟨[(bo†p(t)+bop(t)) (24) ×(ρn,n+100,rvB(t)−ρn+1,n00,rvB(t)) +(brvB†p(t)+brvBp(t))(ρn,n+100,o(t)−ρn+1,n00,o(t))]⟩ = −i∑pκ2p(n+1)ρn00[(nB(ωp)+1)δ(ωp+ω0) (28) +nB(ωp)δ(ωp−ω0)] +i∑pκ2p(n+1)ρn+100[nB(ωp)δ(ωp+ω0) +(nB(ωp)+1)δ(ωp−ω0)],

in which is the Bose-distribution function. With the assumption that the phonon-bath is composed of infinitely many harmonic oscillators having a wide and continuous spectral density, we can make the replacement, in the wide-band limit,

 ∑pκ2p(⋯)⟶∫∞−∞dωpκ2pD(ωp)(⋯), (29)

yielding

 J=−iϖpnB(ω0)(n+1)ρn00+iϖp(nB(ω0)+1)(n+1)ρn+100, (30)

with being constant due to Eq. (3).

Following the same calculational scheme indicated above, we evaluated the other statistical expectation values involved in the EOMs, Eqs. (8c)-(8q). Finally, we obtained the following rate equations in terms of the direct product state representation of the density-matrix for the description of sequential tunneling through a molecular QD, accounting for the unequilibrated phonon effect and its modified tunneling rates, as well as IVM dissipation to the phonon environment:

 ˙ρn00 = ∑m[Γ−1,nmρm11+Γ−2,nmρm22−(Γ+1,nm+Γ+2,nm)ρn00] (31b) −(ϖ+n+ϖ−n)ρn00+ϖ−n+1ρn+100+ϖ+n−1ρn−100, ˙ρn11 = ∑m[Γ+1,mnρm00−(Γ−1,mn+˜Γ+2,nm)ρn11+˜Γ−2,nmρmdd] (31d) −(ϖ+n+ϖ−n)ρn11+ϖ−n+1ρn+111+ϖ+n−1ρn−111, ˙ρn22 = ∑m[Γ+2,mnρm00−(Γ−2,mn+˜Γ+1,nm)ρn22+˜Γ−1,nmρmdd] (31f) −(ϖ+n+ϖ−n)ρn22+ϖ−n+1ρn+122+ϖ+n−1ρn−122, ˙ρndd = ∑m[˜Γ+2,mnρm11+˜Γ+1,mnρm22−(˜Γ−