Time-dependent versus static quantum transport simulations beyond linear response

0.68 MB·English

by  ChiYung Yam

Time-dependent versus static quantum transport simulations beyond linear response
ChiYung Yam,1,2 Xiao Zheng,1 GuanHua Chen,1 Yong Wang,2 Thomas Frauenheim,2 and Thomas A. Niehaus3
1Department of Chemistry, The University of Hong Kong, Pokfulam Road, Hong Kong
2Bremen Center for Computational Materials Science, Am Fallturm 1a, 28359 Bremen, Germany
3University of Regensburg, 93040 Regensburg, Germany∗
(Dated: January 5, 2011)
To explore whether the density-functional theory non-equilibrium Green’s function formalism
(DFT-NEGF)providesarigorousframeworkforquantumtransport,wecarriedouttime-dependent
densityfunctionaltheory(TDDFT)calculationsofthetransientcurrentthroughtworealisticmolec-
ular devices, a carbon chain and a benzenediol molecule inbetween two aluminum electrodes. The
TDDFTsimulationsforthesteadystatecurrentexactlyreproducetheresultsoffullyself-consistent
1
DFT-NEGFcalculationsevenbeyondlinearresponse. Incontrast,sizabledifferencesarefoundwith
1
respect to an equilibrium, non-self-consistent treatment which are related here to differences in the
0
Kohn-Sham and fully interacting susceptibility of the device region. Moreover, earlier analytical
2
conjectures on the equivalence of static and time-dependent approaches in the low bias regime are

n confirmed with high numerical precision.
a
J

4 A first principles description of quantum transport at a)
the atomic scale is usually based on the non-equilibrium

]
i Green’sfunction(NEGF)technique[1],whichreducesto S S
c the traditional Landauer transmission formalism [2] for L R
s
- coherent transport. The device Green’s function is con- Molecular
l
r structedinanad-hocmannerfromtheKohn-Sham(KS) Device
t
m single-particleHamiltonianofgroundstateDensityFunc- L D R
tional Theory (DFT) including the effects of the leads
. b)
t through suitable self-energies [3]. While such methods
a
m provide qualitative results to understand the transport
behavior, there are cases where the theoretical steady
-
d state current differs from experimental results by orders c)
n of magnitude [4]. One reason for the discrepancy is that
o
the evaluated transmission function exhibits resonances
c
[ at the KS single particle energies, which are in general
not physical.
1 FIG. 1. (a) Schematic device setup. (b) Seven-membered
v Recently, several approaches based on time-dependent carbon chain between Al nanowires in (001) direction. Fron-
2 DFT (TDDFT) have been put forward [5, 6], not only tier carbon atoms in hollow position of the Al surface at 1.0
2 to treat genuine dynamical problems like AC transport ˚Adistance. Bondlengths(˚A):C-C=1.32,Al-Al=2.86. (c)
7 or transient effects, but also to remedy the problematic 1,4-benezenediol molecule inbetween Al nanowires.
0
situationinthesteady-statecase. Foropticalabsorption
.
1 spectra, TDDFT excited state energies improve signifi- finally be reached. It is thus argued that TDDFT and
0
cantly on single-particle energy differences that are ob- DFT-NEGFshouldresultintheexactsamesteadystate
1
tained from the static ground state KS equations. Since currentiftheadiabaticapproximationisadoptedforthe
1
: the screening properties of the device depend strongly exchange-correlation functional. We thus encounter a
v
on the optical gap, one would expect also an improved paradox.
i
X description of the current.
A related open question is whether sizable differences
r In fact, dynamical corrections to the conventional between TDDFT and DFT-NEGF may occur for local
a
DFT-NEGFpicturehavebeenpredictedintheregimeof functionals beyond linear response. The interest is here
linear response [7]. Interestingly, these corrections were in the full I-V characteristics of the device including the
showntovanishforapproximatexcfunctionalswhichare bias region around molecular resonances, which is often
localbothintimeandspace. Suchadiabaticfunctionals, used as fingerprint of the contacted molecule in experi-
like the adiabatic local density approximation (ALDA ments. This regime is difficult to access with analytical
a.k.a. TDLDA), are quite common in many implementa- methodsbutamenabletonumericalsimulations. Beyond
tionsofTDDFTforopticalspectraduetotheirnumerical linear response, more than one solution for steady state
simplicity. As shown by Stefanucci et. al. [8] and Zheng current is conceivable, and this was observed for model
et. al. [6], TDLDA leads to the well-known Landauer- systems[9]. Insuchacase,thetransientcurrentmayset-
Bu¨ttiker formula for the steady state current if it can tle into different steady currents depending on how the
2
bias voltage is turned on. Also, it was argued that per-
sistent oscillating currents can exist and a steady state a) 350
mayneverbereached. Thisisbecauseadevicemayhave 300 T=0.1fs T=0.5fs T=2.0fs
bound states with vanishing line widths and correspond- 250 T=0.2fs T=1.0fs
ing infinite life times. The implications of such bound 200
statesforopensystemwereinvestigatedforamodelsys- 150
tem [10], and the current through the model device was 100
found to oscillate for as long as the simulation time. To µA] 50
addresstheseissuesinthisletter,weusearecentlydevel- ent [b) 0
oped TDDFT method for open systems [6] and compare urr
C 40
results with those from DFT-NEGF calculations.
Methodology.— For a two–terminal setup, the entire 30
system consists of a device region (D), and semi-infinite
20
left(L)andright(R)electrodesasshowninFig1a. The
region D, where electron-scattering events take place, is 10
thus an open electronic system.
0
To obtain the reduced electronic dynamics of the de- 0 5 10 15 20
Time [fs]
vice, we focus on the equation of motion (EOM) of σ ,
D
the reduced single-electron density matrix for D, FIG.2. TDLDAtransientcurrentofa)carbonchainandb)
benzenediol for exponential turn-on of the bias voltage with
(cid:88)
iσ˙D =[hD,σD]−i Qα(t). (1) different time constants T and V0 =3V.
α
where h (t) denotes the KS Fock matrix of D, and the effective bias which is smaller than the physical one [7,
D
square bracket on the right–hand side denotes a commu- 8]. In contrast, the current has Landauer form in the
tator. Q isthedissipationtermduetotheαthelectrode ALDA and seems to be equivalent to the static DFT-
α
andbasedontheKeldyshNEGFformalism[1],itcanbe NEGF result. However, the device Green’s function
expressed as,
Gr((cid:15))=[(cid:15)I−h (σ∞)−Σr((cid:15))−Σr((cid:15))]−1, (5)
D D D L R
Q (t)=
α
is constructed from a Hamiltonian that depends on the
(cid:90) t (cid:104) (cid:105)
− dτ Gr(t,τ)Σ<(τ,t)+G<(t,τ)Σa(τ,t) +H.c. steady-state density. The latter is not necessarily the
D α D α
−∞ sameasthedensityobtainedfromaself-consistentDFT-
(2) NEGF treatment or even the equilibrium one. To our
knowledge,ageneraladiabatictheoremforopensystems
whereGr andG< denotethetimedomainretardedand
D D is not available, though adiabaticity has been confirmed
lesser Green’s function in the device region, and Σa and
α for model potentials [11]. Such a theorem would prove
Σ< standfortheadvancedandlesserself-energiesoflead
α the equivalence of TDLDA and LDA-NEGF approaches
α.
from Eq.(5). Hence, we provide a numerical comparison
ThetransientelectriccurrentthroughtheinterfaceS
α of both approaches for some realistic molecular devices
(thecrosssectionseparatingregionD fromtheαthelec-
in this letter.
trode) can be evaluated via [6]:
AsshowninRef.[6],theexpressionforQ (t)inEq.(2)
α
may be considerably simplified in the adiabatic wide–
I (t)=−tr[Q (t)]. (3)
α α bandlimit(AWBL).Itdependsfinallyonlyonthedevice
density matrix σ , the applied bias V(t) and the WBL
Ifthebiaspotentialsapproachaconstantvalueasymp- D
totically (V∞ ≡ V (t → ∞)), a steady state may de- parameters which characterize the lead self-energies [12].
α α
In this way the EOM is closed. At t = 0, the entire
velop. It has been shown that the steady state current
system is assumed to be in thermal equilibrium and the
can then be recast in the form [6, 8]:
initial density matrix σ (t=0) is easily constructed us-
D
(cid:90)
ingequilibriumGreen’sfunction(EGF)techniques. Sub-
I∞ =−I∞ = d(cid:15)[f ((cid:15)−V∞)−f ((cid:15)−V∞)] T((cid:15))
L R L L R R sequently, bias potentials V (t) are applied to the con-
α
T((cid:15))=tr[Gr((cid:15))Γ ((cid:15))Ga((cid:15))Γ ((cid:15))], (4) tacts and uniformly shift the lead energy levels. The
D R D L
device density matrix σ is then propagated according
D
where f are Fermi distribution functions and Γ ((cid:15)) = to Eq.(1), where in each time step the Poisson equation
α α
i[Σr((cid:15))−Σa((cid:15))] describe the lead-molecule coupling. is solved to update the device Hamiltonian.
α α
For non-local functionals, additional XC contributions Results.— We report the calculations of I-V curves for
appear in the lead Fermi functions and give rise to an two molecular systems, a carbon chain and a benzene-
3
state current is the same in all cases. This was not
obvious, as the history of the applied bias voltage is
a) 350 00..45 IT∞DLDA ILNDEAGF ILEDGAF included in the formalism, although the memory effect
µent [A] 223050000 000...123µCurrent [A] oA((cid:39)fLtD1h0Aef.se)Txachnhedaenosgtcaec-buclroissrhirnemleaevtneitornyofcftauhsneec.stWtioenaedadlyoisnstoaattbeosbiessneutrlvtirneapftaehsret-
Curr 150 0 1 Vo l2tage [ m3V] 4 5 0 sistentcurrentoscillations[10],asthetermQα inEq.(2)
100
ensures proper dissipation in the leads.
50
TDLDA simulations were then performed for several
0
0 1 2 3 4 5 bias values in the linear response, low voltage, regime.
b) 1 Voltage [V] Theasymptoticvaluesofthetransientcurrents(I∞ )
gth 0.8 FKuolhl nre-Ssphoanmse were extracted in each case and compared to statTicDDLDFAT
pole Stren[arb. unit] 000 ...0246 ctiaolncufloartitohneslienatdhseinLDthAeuinsisnegtstohfeFsaigm.e3(Wa)BaLndapFpirgo.x4im(aa)-.
Di 0 2 4 6 8 10 TwokindsofstaticDFTresultsarereported. Inthefirst
Energy [eV]
case,thedeviceGreen’sfunctionisdeterminedfullyself-
FIG. 3. (a) I-V curves of seven-membered carbon chain be-
consistently using the non-equilibrium density (INEGF).
tween Al nanowires in (001) direction at high bias and low LDA
bias (inset) regime. (b) TDLDA optical spectrum of isolated In the second case, the equilibrium density is employed
and hydrogen capped carbon chain (C H ) (see text). (IEGF) and the Poisson equation is solved only once,
7 2 LDA
which yields a linear potential drop across the molecule.
diolmoleculesandwichedbetweenaluminumleadsinthe
For both studied systems all approaches yield identical
(001) direction of bulk Al. The structures are shown
results. This provides numerical evidence for the earlier
in Fig. 1(b) and 1(c) respectively. Simulations for an
analytical arguments [7] on the equivalence of static and
Al wire were also performed with similar findings as re-
dynamic local DFT approaches in linear response.
ported below. For each of the molecular systems, we in-
In another set of simulations, the covered bias range
cludeexplicitly18Alatomsofeachcontactinthedevice
wasextendedto5V.Finitebiassimulationsallowforthe
region to ensure smooth potentials. The minimal Gaus-
coverage of the resonant tunneling regime in which the
sianbasissetSTO-3Gisadoptedinthecalculations. For
biaswindowcontainsoneormoreofthemolecularlevels.
the propagation, the fourth-order Runge-Kutta method
Inspectionofthetransmissionfunction(notshownhere)
isemployedwithatimestepof0.02fstointegrateEq.(1)
shows that transport occurs predominantly through the
in the time domain. The ALDA is adopted for exchange
unoccupied levels in both systems. For benzenediol, the
and correlation. The bias voltage is turned on expo-
conducting molecular orbitals (MO) closest to the Fermi
nentially at t = 0 with a time constant T, according
energy (E ) are located 0.92 eV above and 2.58 eV be-
F
to V(t)=V (1−e−t/T).
0 low E , respectively. The carbon chain exhibits metallic
F
behavior with a partially filled lowest unoccupied MO
a) 100 25 I∞ INEGF IEGF derived level [13].
20A] TDLDA LDA LDA ThefullI-Vcharacteristicsofbothsystemsaresumma-
µnt [A] 6800 51105Current [n arinzdedLinDAF-igN.E3G(aF) a(nILdNDEFAGiFg). 4c(uar)r.enTthetrTacDeLsDaAre(vIT∞irDtuLDalAly)
urre 40 0 1 2 3 4 5 0 identical over the whole bias range including the region
C Voltage [mV] of resonant transmission. This level of agreement sup-
20 ports not only the validity of the Landauer picture in
thecontextofTDLDA(Eq.(4)), butalsoshowsthatthe
0
b) 0 1 2Voltage [V] 3 4 5 TDLDA and LDA-NEGF densities are in fact identical
1 forrealisticdevicesevenforstrongbiasvoltage. Another
h Full response
pole Strengt[arb. unit] 0000 ....02468 Kohn-Sham rcNeuEarsrGoeFnnctetoroesfatmtumdueyltnittp.hleeAshosilgudhtiisbocniuassssienrdetghibmeyseeSli´afs-nctcohhneeszipsoteestnsiatblLl.eD[oA9c]--,
Di 2 3 4 5 6 7 8 9 10
Energy [eV] a dynamical approach would single out the most stable
FIG.4. (a)I-Vcurvesof1,4-benzenediolinbetweenAlleads steady-state solution in such a situation. Attempts to
at high bias and low bias (inset) regime. (b) TDLDA optical find such multiple solutions were unsuccessful for both
spectrum of isolated 1,4-benzenediol (C6O2H6). devices. We also performed simulations with increased
molecule-leaddistanceinordertocoveralargerparame-
Fig.2depictsthetransientcurrentsforthetwomolec- terspaceofcouplingmatrixelements. Thetransientcur-
ular systems at different values of T with V = 3 V. Al- rents show oscillations with slower decay to the steady-
0
though the initial current traces differ, the final steady- stateinthesecases,butagain,theTDLDAresultsmatch
4
the corresponding LDA-NEGF simulations exactly. We thank the German Science Foundation (DFG,
Comparing now the TDLDA or LDA-NEGF results SPP 1243), Hong Kong Research Grant Council
with LDA-EGF values and focussing first on the bias re- (HKU700909P, 700808P, 701307P), Hong Kong Univer-
gion below 2 V, we find significant differences which are sity Grant Council (AoE/P-04/08) and National Science
more pronounced in the carbon chain. The LDA-EGF Foundation of China (NSFC 20828003) for support.
current is initially higher, which can be traced back to a
movement of the LUMO resonance to higher energies in
theself-consistentapproaches. Infact,TDDFTexplicitly
includesthepotentialchangeduetotheinduceddensity,
∗ Corresponding author: [e-mail],
similar to the random-phase-approximation. In the con-
[e-mail]
text of harmonic perturbations in the optical range, this
[1] L. P. Kadanoff and G. Baym, Quantum statistical me-
leads to the well known shift of TDDFT excited states
chanics: Green’s function methods in equilibrium and
away from simple KS single-particle energy differences. nonequilibrium problems (Addison-Wesley, 1989).
In order to quantify this effect, we computed the optical [2] R. Landauer, IBM J. Res. Dev. 1, 233 (1957);
spectrum of the isolated benzenediol and carbon chain M. Bu¨ttiker, Physical Review Letters 57, 1761 (1986).
using a linear response TDDFT method [14]. Fig. 3(b) [3] J. C. Cuevas and E. Scheer, Molecular Electronics: An
IntroductiontoTheoryandExperiment(WorldScientific,
and 4(b) depict the results using the full response and a
2010).
non-self-consistenttreatmentemployingthegroundstate
[4] S.M.LindsayandM.A.Ratner,AdvancedMaterials19,
density, i.e. the KS response. Although in both systems
23 (2007).
the KS and TDDFT spectra differ strongly, the devi- [5] J. K. Tomfohr and O. F. Sankey, Phys. stat. sol. (b)
ations are higher in the carbon chain with a shift of 226, 115 (2001); R. Baer, T. Seideman, S. Ilani,
roughly 3 eV compared to a shift of 0.5 eV for benzene- and D. Neuhauser, J. Chem. Phys. 120, 3387 (2004);
diol. This finding is in line with the transport results K. Burke, R. Car, and R. Gebauer, Phys. Rev. Lett.
94, 146803 (2005); S. Kurth, G. Stefanucci, C. O. Alm-
for V < 2 V. For larger values of the bias, the observed
bladh, A. Rubio, and E. K. U. Gross, Phys. Rev. B 72,
current traces are more difficult to interpret since the
35308(2005); J.Yuen-Zhou,C.Rodr´ıguez-Rosario, and
transmissionprofileisalteredinanon-linearfashionwith
A. Aspuru-Guzik, Phys. Chem. Chem. Phys. 11, 4509
respect to the equilibrium situation. (2009).
Summary.—In this letter, we performed TDLDA [6] X.Zheng,F.Wang,C.Y.Yam,Y.Mo, andG.H.Chen,
transportsimulationsbasedonatimepropagationofthe Phys. Rev. B 75, 195127 (2007).
KS density matrix and compared these to conventional [7] F. Evers, F. Weigend, and M. Koentopp, Phys. Rev.
B 69, 235411 (2004); N. Sai, M. Zwolak, G. Vignale,
static DFT results in the Landauer picture. The studied
and M. Di Ventra, Phys. Rev. Lett. 94, 186810 (2005);
systems comprise the ubiquitous benzenediol molecule
G. Vignale and M. Di Ventra, Phys. Rev. B 79, 14201
withsmallequilibriumconductanceandanatomicchain
(2009).
with metallic conduction properties. In the regime of [8] G. Stefanucci, S. Kurth, E. K. U. Gross, and A. Rubio,
linear response, the analytically predicted equivalence of Theo. and Comp. Chem. 17, 247 (2007).
static and time-dependent DFT formalisms for local xc [9] C.G.Sa´nchez,M.Stamenova,S.Sanvito,D.R.Bowler,
potentials could be confirmed with high numerical preci- A.P.Horsfield, andT.N.Todorov,J.Chem.Phys.124,
214708 (2006).
sion. Beyondlinearresponse,TDLDAreproducesthere-
[10] E.Khosravi,S.Kurth,G.Stefanucci, andE.K.U.Gross,
sult of LDA-NEGF calculations and deviates from static
Appl. Phys. A 93, 355 (2008).
LDA only when the equilibrium density is employed in
[11] G.ScheitlerandM.Kleber,Zeit.f.Phys.D9,267(1988).
the latter. This parallels the energy shifts seen in go- [12] Alternatively (but not employed in this study), a hier-
ing from the KS response to the full coupled response in archical equation of motion approach [17] was recently
DFTabsorptionspectra. Althoughthetransientcurrent developed,forwhichthehierarchyterminatesexactlyat
dependsonthespecificsofappliedbiasvoltage,aunique the second tier without any approximation.
[13] B.Larade,J.Taylor,H.Mehrez, andH.Guo,Phys.Rev.
steadystateisreachedforthesameeventualbiasvoltage.
B 64, 75420 (2001).
More advanced xc functionals beyond the ALDA are
[14] C. Y. Yam, S. Yokojima, and G. H. Chen, Physical Re-
expectedtoleadtoanimprovementinvariousways. For
view B 68, 153105 (2003).
example, non-local functionals provide sizable dynami- [15] C.Toher,A.Filippetti,S.Sanvito, andK.Burke,Phys.
cal corrections to the Landauer formula [7]. Moreover, Rev.Lett.95,146402(2005); S.H.Ke,H.U.Baranger,
functionals which exhibit a proper derivative disconti- and W. Yang, J. Chem. Phys. 126, 201102 (2007).
nuity refine not only the energetical position of molecu- [16] C.Y.Yam,Y.Mo,F.Wang,X.Li,G.H.Chen,X.Zheng,
Y.Matsuda,J.Tahir-Kheli, andA.G.WilliamIII,Nan-
lar levels, but also the molecule-lead coupling [15]. De-
otech. 19, 495203 (2008).
spiteoftheseshort-comings,TDLDAisstillusefulinthe
[17] X.Zheng,G.H.Chen,Y.Mo,S.K.Koo,H.Tian,C.Y.
qualitative description of transport problems which re-
Yam, andY.J.Yan,J.Chem.Phys.133,114101(2010).
quire a dynamical treatment, like the transient behavior
of molecular devices or photo-switches [16].