Table Of ContentIterative summation of path integrals for nonequilibrium molecular quantum transport
R. Hützen,1 S. Weiss,1,2 M. Thorwart,2 and R. Egger1
1Institut für Theoretische Physik, Heinrich-Heine-Universität, D-40225 Düsseldorf, Germany
2I. Institut für Theoretische Physik, Universität Hamburg, D-20355 Hamburg, Germany
(Dated: February 1, 2012)
We formulate and apply a nonperturbative numerical approach to the nonequilibrium current,
I(V), through a voltage-biased molecular conductor. We focus on a single electronic level coupled
to an unequilibrated vibration mode (Anderson-Holstein model), which can be mapped to an ef-
fective three-state problem. Performing an iterative summation of real-time path integral (ISPI)
expressions, we accurately reproduce known analytical results in three different limits. We then
2
studythecrossoverregimebetweenthoselimitsandshowthattheFranck-Condonblockadepersists
1
inthequantum-coherentlow-temperaturelimit,withanonequilibriumsmearingofstepfeaturesin
0
the IV curve.
2
n PACSnumbers: 73.63.-b,72.10.-d,02.70.-c
a
J
1 Introduction.—Understanding quantum transport in is also based on a Keldysh path-integral formulation but
3 nanoscaleelectronicsystemswithvibrationalormechan- does not involve stochastic sampling schemes and thus
ical (“phonon”) degrees of freedom is of topical interest remains free from any sign problem. To that end, we ex-
]
in several areas of physics, including molecular electron- tend the “iterative summation of path integrals” (ISPI)
l
l
a ics [1, 2], inelastic tunneling spectroscopy [3], nanoelec- technique[19]totheAHmodel. Technicalaspectsofthe
h tromechanical systems [4], break junctions [5], and sus- present approach, in particular our mapping to an effec-
- pended semiconductor or carbon-based nanostructures tive three-state system via the ”spin-1 Hirsch-Fye trans-
s
e [6–9]. The electron-phonon interaction allows to observe formation” in Eq. (2) below, should also be of interest
m
a rich variety of intriguing phenomena, such as negative to QMC schemes [14]. In essence, the ISPI method ex-
. differential conductance, the Franck-Condon blockade of ploits that time correlations of the auxiliary three-state
t
a transport, rectification, vibrational sidebands or steplike Keldyshvariable,whichariseafterfunctionalintegration
m
features in the current-voltage (IV) characteristics, and over the phonon and the (dot and lead) fermion degrees
- current-induced heating or cooling. Already the sim- of freedom, can be truncated beyond a certain memory
d
plest nonequilibrium “Anderson-Holstein” (AH) model, time τ when either T or V is finite. Together with a
n m
o where the nanostructure corresponds to just one spin- convergence scheme designed to eliminate systematic er-
c less electronic level coupled to a single oscillator mode, rorsduetothefinitenessofτ ,suchcalculationsallowto
m
[ captures much of this richness [10, 11]; for a review, see obtain numerically exact results. The ISPI method has
1 Ref. [3]. Analytical approaches allow to understand the already been successfully applied to the spinful Ander-
v AH model in various corners of parameter space, but son model [19–21], where instead of the phonon a local
6 no controlled approximation, let alone exact solution, charging interaction is present. While we focus on the
6
connecting these corners seems in reach. One may ex- simplest version of the AH model with a single unequi-
4
pect that a unified picture is available from numerics. librated [22] phonon mode here, the conceptual gener-
6
. However, numerical renormalization group [12] or quan- alization to include Coulomb interactions, more phonon
1
tum Monte Carlo (QMC) calculations [13, 14] are usu- modes,orseveraldotlevelsisstraightforward. Webench-
0
ally restricted to equilibrium. For the nonequilibrium markourISPIcodeagainstthreedifferentanalyticalap-
2
1 AH model, Han [15] employed an imaginary-time QMC proaches and then study the crossover between the re-
: approach followed by a double analytical continuation spective regimes.
v
scheme; unfortunately, the latter step is plagued by in- AH Model.—We consider the AH Hamiltonian, H =
i
X
stabilities [16]. A promising avenue for the AH model H + H + H , describing a molecular level with tun-
m t l
r has recently been suggested by real-time path-integral nel coupling (H ) to metallic source and drain contacts.
a t
QMC simulations [17, 18], where one directly computes Taking a single spinless dot level (fermion annihilation
the time-dependent current. Such calculations have to operator d) with energy (cid:15) and a boson mode (annihila-
deal with the infamous dynamical sign problem at long tion operator b) of frequency Ω, the isolated molecule
times,butinseveralparameterregions,especiallywhena Hamiltonian is [3] (we use units with (cid:126)=k =1)
B
secondaryphononbathispresent,thestationarysteady-
state regime can be reached. Hm =Ω b†b+(cid:2)(cid:15)+λ(b+b†)(cid:3)nd (1)
Inthiswork,weformulateandapplyanalternativenu- withn =d†dandtheelectron-phononcouplingstrength
d
merical approach, which in practice is useful unless both λ. The lead Hamiltonian H is taken in the standard
l
the temperature T and the bias voltage V are small. It wide-band approximation [23], with relaxation processes
2
assumed fast enough to have Fermi functions with tem- Keldysh path integral are then effectively noninteracting
peratureT aselectronicdistributions; theirchemicalpo- butcoupletothetime-dependentauxiliaryspinvariable.
tential difference defines the bias voltage V. The tun- Hencethosefieldscanbeintegratedoutanalyticallyand
nel coupling then introduces the hybridization energy the time-dependent current, I(t ), follows from a path-
j
scales Γ and Γ for the left/right lead [23]. For sim- summation formula for the generating functional,
L R
plicity, we focus on the symmetric case in what follows,
(cid:88)
Γ = Γ = Γ/2. As observable of main interest, we Z = detD[{s}], (3)
L R
study the steady-state current I through the molecule. {sη=0,±1}
Analytical approaches.—Before turning to a descrip-
where the matrix D (in time and Keldysh space) de-
tion of the ISPI scheme, let us briefly summarize the ηη(cid:48)
pends on the complete spin path {s}. Specifically, we
analytical approaches to the AH model that we employ
obtain D =−iB(G−1−Σ), where G−1 is the discretized
to benchmark our method. (i) For λ/Γ (cid:28) 1, perturba- d d
inverse Green’s function of the dot as in Refs. [19, 31]
tion theory in the electron-phonon coupling applies and
but with the modified spin-dependent matrix elements
yields a closed IV expression for arbitrary values of all (cid:2)−iG−1(cid:3) = −s . The self-energy matrix Σ de-
other parameters [24]. We note that the solution of the d η+1,η η
scribesthetraced-outleads. WefindΣ (cid:54)=0onlywhen
AH model with a very broad dot level [25, 26] corre- ηη(cid:48)
s = ±1, where it coincides with the usual (wide-band
spondstothissmall-λregime. (ii)Forhightemperatures, η
limit) expression [23]. Finally, the diagonal matrix B
T (cid:29) Γ, a description in terms of a rate equation is pos-
(quoted here for (cid:15)=0) with
sible [23]. We here use the simplest sequential tunneling
rveesrpsioonndwinigthregsouldltesnmrualtechratthesos[e27o].f Fpoerrtsumrbaalltiλo,nththeecoorry-, Bηη =Asηe−λ2δt2(cid:80)η(cid:48)σσ(cid:48)[iGph]η,η(cid:48)+1|sηsη(cid:48)| (4)
while in the opposite strong-coupling limit, the Franck- encapsulates all phonon effects, where G is the dis-
ph
Condon blockade occurs and implies a drastic current cretized phonon Green’s function, see Ref. [31], and we
suppression at low bias voltage [8, 28]. (iii) For small os- used the notation A0 = 1 and A±1 = ±(1/2)e−λ2δt2/2.
cillator frequency, Ω (cid:28) min(Γ,eV), the nonequilibrium By including a source term in D, it is straightforward
Born-Oppenheimer (NEBO) approximation is controlled to numerically extract the time-dependent current I(t )
j
and allows to obtain I from a Langevin equation for the [19]. Forsufficientlylongtimest ,I(t )reachesaplateau
j j
oscillator [29, 30]. For small λ, this approach is also con- yielding the steady-state current of interest.
sistentwithperturbativetheory,whileforhighT,NEBO ISPI implementation.—Starting from the formally ex-
and rate equation results are found to agree. act path-summation formula [Eq. (3)], the ISPI algo-
Keldysh path-summation formulation.—We now start rithm can now be adapted from Ref. [19]: For finite T
from the textbook time-discretized coherent-state repre- orV,matrixentriesD involvinglargetimedifferences
ηη(cid:48)
sentation of the Keldysh generating functional [31, 32]. |t−t(cid:48)| are exponentially small. We put these to zero be-
The short-time propagator on the forward/backward yond a memory time τ ≡ Kδ , where K denotes the
m t
branch of the Keldysh contour, e∓iδtH, where δt denotes number of time slices kept in the memory. The numeri-
thediscretetimestep,thenallowsforaTrotterbreakup, cal computation of the memory-truncated path summa-
e∓iδtH = e∓iδtH1e∓iδt(H−H1), with the systematic er- tion in Eq. (3) is then possible in an iterative way [19]
ror in observables scaling ∼ δt2. It is useful to choose without additional approximations and, for given δt and
H1 = Hm−Ωb†b, see Eq. (1), where the auxiliary rela- K, yields the steady-state current I(δt,K). While this
tion formulation is exact for δ → 0 and sufficiently large
t
e∓iδtH1 =1−nd+nde−λ2δt2/2e∓iδt(cid:15)e∓iδtλb†e∓iδtλb (2) K (long memory time), K cannot be chosen arbitrarily
large in practice and an extrapolation scheme is neces-
allows to effectively decouple the electron-phonon inter- sary[19]. Convergenceoftheextrapolationrequiressuffi-
action in terms of a three-state variable s = 0,±1 de- ciently high T or V, for otherwise the necessary memory
η
fined at each (discretized) time step t along the for- times are exceedingly long. For the results below, we
j
ward/backward (σ = ±) part of the Keldysh contour, used K ≤ 4 and 0.3 ≤ Γδ ≤ 0.35. The shown current
t
where η = (t ,σ). Below, we also use the notation follows by averaging over the δ -window, with error bars
j t
η±1 = (t ,σ) with periodic boundary conditions on indicating the mean variance. Additional ISPI runs for
j±1
the Keldysh contour. It is crucial for the construction of 0.18 ≤ Γδ ≤ 0.22 and 0.3 ≤ Γδ ≤ 0.4 were consistent
t t
thecoherent-statefunctionalintegralthatEq.(2)isnor- withtheseresults, andweconcludethatsmallerrorbars
mal ordered. The “spin” variable s picks up the three indicate that convergence has been reached. For typical
η
terms in Eq. (2) and acts like a Hubbard-Stratonovich parameters and K = 4, our ISPI code yielding I(δ ,K)
t
auxiliary field, similar to the Ising field employed in the runs for ≈11 CPU hours on a 2.93 GHz Xeon processor.
Hirsch-Fye formulation of the Anderson model [19, 33]. Benchmark checks.—Nextweshowthatthenumerical
The bosonic (phonon) scalar field and the fermionic (dot ISPI results are consistent with analytical theory for the
and lead electrons) Grassmann fields appearing in the IV curvesinallthreeparameterlimitsmentionedabove.
3
1.20.4 T = 3Γ IpSePrtI., tThe=oΓry 0.6 0.3 INrSaEtPeBI e,O qT.=3Γ INSEPBI,O T = Γ
1 rate eq. 0.2
0.2
0.1
)0.8 )
Γ/h 00 0.5 1 1.5 2 Γ/h0.4 00 0.5 1 1.5 2
(e0.6 (e
I / I / 0.8 T = Γ/3
1.5 0.6
0.4 T = Γ/3 0.2
1 0.4
0.2 0.5 0.2
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
0 0
0 0.5 1 1.5 2 0 0.5 1 1.5 2
eV / Γ eV / Γ
Figure 1. (Color online) Current I (in units of eΓ/h) vs bias Figure3. (Coloronline)SameasFig.1butforΩ=0.5Γand
voltage V (in units of Γ/e) for λ = 0.5Γ, Ω = Γ, (cid:15) = 0, λ = Γ. The main panel is for T = Γ and compares the ISPI
and T = Γ. The ISPI data are depicted as filled red circles, results to NEBO predictions. The insets are for T =3Γ and
wherethedottedredcurveisaguidetotheeyeonlyandthe T = Γ/3, respectively, where also the rate equation results
error bars are explained in the main text. We also show the areshown. NoticethatincontrasttoISPI,therateequation
results of perturbation theory in λ (dashed black curve) and predicts an unphysical current blockade for T =Γ/3.
of the rate equation (solid blue curve). The upper (lower)
inset shows the corresponding result for T =3Γ (T =Γ/3).
memory time (τ ) demands. The effect of changing the
m
phononfrequencyΩisillustratedinFig.2,takingT =Γ
but otherwise identical parameters. Again perturbation
theory is well reproduced. Next, Fig. 3 shows ISPI re-
1.2 1 Ω = 2Γ ISPI, Ω = Γ/2 sults for a slow phonon mode, Ω = Γ/2, with stronger
pert. theory electron-phonon coupling, λ = Γ. In that case, pertur-
1 rate eq.
bation theory in λ is not reliable and the rate equation
0.5
is only accurate at the highest temperature (T = 3Γ)
)0.8
h
Γ/ studied, cf. the upper left inset of Fig. 3. However, we
(e0.6 00 0.5 1 1.5 2 observe from Fig. 3 that for such a slow phonon mode,
I / NEBO provides a good approximation for all tempera-
0.4
tures and/or voltages of interest. We conclude that the
ISPI technique is capable of accurately describing three
0.2
different analytically tractable parameter regimes.
0 Franck-Condon (FC) blockade.—Next we address the
0 0.5 1 1.5 2
eV / Γ limit of strong electron-phonon coupling λ, where the
rate equation approach yields a FC blockade of the cur-
rentforlowbiasandT (cid:29)Γ[28]. Sufficientlylargeλcan
Figure 2. (Color online) Same as Fig. 1 but for Ω = 0.5Γ be realized experimentally, and the FC blockade has in-
(main panel) and Ω=2Γ (inset), both for T =Γ. deedbeenobservedinsuspendedcarbonnanotubequan-
tum dots [8]. For an unequilibrated phonon mode with
intermediate-to-large λ, understanding the FC blockade
For clarity, we focus on a resonant level with (cid:15)=0 here. inthequantum-coherentlow-temperatureregime,T <Γ,
Let us then start with the case of weak electron-phonon is an open theoretical problem. Here multiple phonon
coupling, λ = 0.5Γ. Figure 1 compares our ISPI data excitation and deexcitation effects imply a complicated
for Ω = Γ to the respective results of perturbation the- (unknown)nonequilibriumphonondistributionfunction,
ory in λ and of the rate equation. As expected, for this andtheone-steptunnelinginterpretationintermsofFC
parameter choice, perturbation theory essentially repro- matrixelementsbetweenshiftedoscillatorparabolas[28]
duces the ISPI data. The rate equation is quite accurate is not applicable anymore. We here study this question
for high temperatures, but quantitative agreement with usingISPIsimulations,whichautomaticallytakeintoac-
ISPI was obtained only for T (cid:38)10Γ. Note that the ISPI count quantum coherence effects.
error bars increase when lowering T due to the growing In Fig. 4, the crossover from weak to strong electron-
4
[2] N.A. Zimbovskaya and M.R. Pederson, Phys. Rep. 509,
1 (2011).
[3] M. Galperin, M.A. Ratner, and A. Nitzan, J. Phys.:
λ = 0.5Γ Τ=0.2Γ
Cond. Matt. 19, 103201 (2007).
T = Γ λλ == 12..55ΓΓ [4] H.G.Craighead,Science290,1532(2000);M.L.Roukes,
1
3 λ = 3Γ Phys. World 14, 25 (2001).
λ = 4Γ [5] O.Tal,M.Krieger,B.Leerink,andJ.M.vanRuitenbeek,
h) Phys. Rev. Lett. 100, 196804 (2008).
Γ/I / (e2 00 2 4 6 8 [6] GMSc.eiAeenr.wcSeatle3de2tl,e5,L, .A1P1..0KK3.o(Hu20wü0tet9ne)hl.,ovBe.n,WaintdkaHm.Sp.,JM.v.anPodoetr,ZHa.nBt.,
1 [7] B. Lassagne, Y. Tarakanov, J. Kinaret, D. Garcia-
Sanchez, and A. Bachtold, Science 325, 1107 (2009).
[8] R. Leturcq, C. Stampfer, K. Inderbitzin, L. Durrer, C.
Hierold, E. Mariani, M.G. Schultz, F. von Oppen, and
0
0 2 4 6 8 K. Ensslin, Nature Phys. 5, 327 (2009).
eV / Γ
[9] A. K. Hüttel, B. Witkamp, M. Leijnse, M.R. Wegewijs,
and H.S.J. van der Zant, Phys. Rev. Lett. 102, 225501
(2009).
Figure 4. (Color online) ISPI data for the IV curves from [10] A. Mitra, I. Aleiner, and A.J. Millis, Phys. Rev. B 69,
weak(λ=0.5Γ)tostrong(λ=4Γ)electron-phononcoupling, 245302 (2004).
with Ω = 2Γ. The main panel is for T = 0.2Γ, the inset for [11] H. Wang and M. Thoss, J. Chem. Phys. 131, 024114
T = Γ. We used a dense voltage grid yielding smooth IV (2009).
curves. Errorbarsarenotshownbutremainsmall,cp.Fig.1. [12] P.S. Cornaglia, H. Ness, and D.R. Grempel, Phys. Rev.
Lett. 93, 147201 (2004).
[13] L. Arrachea and M.J. Rozenberg, Phys. Rev. B 72,
041301(R) (2005).
phonon coupling λ is considered. The inset shows IV
[14] E.Gull,A.J.Millis,A.I.Lichtenstein,A.N.Rubtsov,M.
curvesforT =Γ,whereweobserveacurrentblockadefor Troyer,andP.Werner,Rev.Mod.Phys.83,349(2011).
low voltages once λ (cid:38) 2Γ. The blockade becomes more [15] J.E. Han, Phys. Rev. B 81, 113106 (2010).
pronounced when increasing λ and is lifted for voltages [16] A. Dirks, P. Werner, M. Jarrell, and T. Pruschke, Phys.
abovethepolaronenergyλ2/Ω[28]. Remarkably,theFC Rev. E 82, 026701 (2010).
[17] L. Mühlbacher and E. Rabani, Phys. Rev. Lett. 100,
blockadepersistsandbecomesevensharperasoneenters
176403 (2008).
the quantum-coherent regime (here, T = 0.2Γ), despite
[18] M. Schiró and M. Fabrizio, Phys. Rev. B 79, 153302
ofthebreakdownofthesequentialtunnelingpicture. We
(2009).
also observe a nonequilibrium smearing of phonon step [19] S. Weiss, J. Eckel, M. Thorwart, and R. Egger, Phys.
features in the IV curves in Fig. 4, cf. also Refs. [8, 28]. Rev. B 77, 195316 (2008); ibid. 79, 249901(E) (2009).
[20] J.Eckel,F.Heidrich-Meisner,S.G.Jakobs,M.Thorwart,
Conclusions.—We have extended the iterative simula-
M. Pletyukhov, and R. Egger, New J. Phys. 12, 043042
tion of path integrals (ISPI) technique to the Anderson-
(2010).
Holstein model, which is the simplest nonequilibrium
[21] D. Segal, A.J. Millis, and D.R. Reichman, Phys. Rev. B
model for quantum dots or molecules with an intrinsic 82, 205323 (2010).
bosonic(phonon)mode. Ourformulationexploitsamap- [22] For an equilibrated mode, i.e., in the presence of a sec-
ping to an effective three-state system and reproduces ondary phonon bath, both ISPI and other numerical
three analytical theories valid in different parameter re- schemes [11, 17, 18] simplify substantially.
[23] Yu.V. Nazarov and Ya.M. Blanter, Quantum Transport
gions. ThisextensionoftheISPIapproachthencaptures
(Cambridge University Press, 2010).
thefullcrossoverbetweenthoselimitsunlessbothT and
[24] R. Egger and A.O. Gogolin, Phys. Rev. B 77, 113405
V are very small. For strong electron-phonon coupling
(2008).
and an unequilibrated phonon mode, we find that the [25] B. Dorá and A. Halbritter, Phys. Rev. B 80, 155402
Franck-Condonblockadebecomesevenmorepronounced (2009).
as one enters the quantum coherent regime. [26] Y.Vinkler,A.Schiller,andN.Andrei,Phys.Rev.B85,
035411 (2012).
WethankA.Zazunovforhelpfuldiscussions. Financial
[27] S. Braig and K. Flensberg, Phys. Rev. B 68, 205324
supportbytheDFG(SPP1243andSFB668)andbythe
(2003).
ZIM (Düsseldorf) is acknowledged.
[28] J. Koch and F. von Oppen, Phys. Rev. Lett. 94, 206804
(2005).
[29] F. Pistolesi, Ya.M. Blanter, and I. Martin, Phys. Rev. B
78, 085127 (2008).
[30] N.Bode,S.ViolaKusminskiy,R.Egger,andF.vonOp-
[1] J.C. Cuevas and E. Scheer, Molecular electronics: An pen, Phys. Rev. Lett. 107, 036804 (2011).
IntroductiontoTheoryandExperiment(WorldScientific, [31] A. Kamenev and A. Levchenko, Adv. Phys. 58, 197
Singapore, 2010). (2009).
5
[32] A. Altland and B.D. Simons, Condensed Matter Field [33] J.E. Hirsch and R.M. Fye, Phys. Rev. Lett. 56, 2521
Theory, 2nd ed. (Cambridge University Press, 2010). (1986).
