Table Of Content

A time-dependent approach to electron pumping in open quantum systems G. Stefanucci,1,2,∗ S. Kurth,1,2 A. Rubio,2,3 and E. K. U. Gross1,2 1Institut fu¨r Theoretische Physik, Freie Universität Berlin, Arnimallee 14, D-14195 Berlin, Germany 2European Theoretical Spectroscopy Facility (ETSF) 3Departamento de F´ısica de Materiales, Facultad de Ciencias Qu´ımicas, 7 UPV/EHU, Unidad de Materiales Centro Mixto CSIC-UPV/EHU and 0 Donostia International Physics Center (DIPC), San Sebastián, Spain 0 (Dated: February 4, 2008) 2 n We propose a time-dependentapproach to investigate the motion of electrons in quantum pump a deviceconfigurations. Theoccupiedone-particlestatesarepropagatedinrealtimeandusedtocal- J culatethelocal electron densityandcurrent. Anadvantageofthepresentcomputationalschemeis 2 thatthesamecomputationaleffortisrequiredtosimulatemonochromatic, polychromaticandnon- 1 periodicdrivings. Furthermore,initialstatedependenceandhistoryeffectsarenaturallyaccounted for. This approach can also be embedded in the framework of time-dependent density functional ] theory to include electron-electron interactions. In the special case of periodic drivings we com- ll binethe Floquet theory with nonequilibrium Green’s functions and obtain a general expression for a the pumped current in terms of inelastic transmission probabilities. This latter result is used for h benchmarkingour propagation scheme in thelong-time limit. Finally, we discuss thelimitations of - s Floquet-based schemes and suggest our approach as a possible way to go beyond them. e m PACSnumbers: 05.60.Gg,72.10.-d,73.23.-b,73.63.-b . t a I. INTRODUCTION speaking,Floquet-basedapproachesprovidea very pow- m erful tool to calculate average quantities of periodically - driven systems. However, going beyond the monochro- d The continuous progress in manipulating single n moleculeschemicallyboundtomacroscopicreservoirshas matic case quite quickly becomes computationally de- o led to the emerging field of molecular electronics.1 Be- manding. Furthermore, such approaches are not appli- c cable to the study of transient effects and non-periodic sides the widely studied stationary case, today experi- [ mental techniques enable the study of time-dependent phenomena. 1 phenomena in open quantum systems, like photon- v assistedtransportandelectronpumpingthroughrealistic In this work we propose a time-dependent approach 9 or artificial molecules. suited to study the effects of an electric field, like a gate 7 2 An electron pump is an electronic device generating a voltage or a laser field, on the electron dynamics of a 1 netcurrentbetweentwounbiasedelectrodes. Pumpingis nanoscale junction. Our approach allows for calculat- 0 achieved by applying a periodic gate voltage depending ing the full time dependence (including the transientbe- 7 on two or more parameters. Electron pumps have been havior) of observables like the local density or current, 0 realized experimentally, e.g., for an open semiconductor and the same computational effort is required for both / t quantumdot2 drivenbytwoharmonicgatevoltageswith monochromaticandpolychromaticdrivingsaswellasfor a m a phase shift, and for a open nanotube3 driven by an nonperiodic perturbations. electrostatic potential wave. - d Intheliteraturedifferenttechniqueshavebeenusedto The paperis organizedasfollows. InSectionII wede- n discuss electron pumping theoretically. For slowly vary- scribethesystemconsistingoftwomacroscopicreservoirs o ing electric field the device remains in equilibrium and connected to a central device. We combine the Floquet c the pumping process is adiabatic.2,4 Brouwer5 has sug- theorywiththeKeldyshformalismtostudythelong-time : v gestedascatteringapproachtodescribeadiabaticpump- behavior of the device, and we generalize the formula Xi ings, but his treatment is limited to periodic potentials for the average current by Camalet et al..7 Some gen- depending on only two parameters. The generalization eralfeatures of a Floquet-based algorithmare discussed. r a to arbitrary periodic potentials has been put forward by To overcome the limitations of the Floquet theory we Zhouetal.6 whousedtheKeldyshtechniquetocalculate describe a real-time approach based on the propagation the net charge transported across the device per period. of the occupied single-particle states. Full implementa- A natural way to go beyond the adiabatic case is to tion details are given for one-dimensional electrodes and apply Floquet theory. Within an equation-of-motion arbitrary device geometries. The performance of the al- approach Camalet et al.7 have found a general expres- gorithmisillustratedinSectionIIIwherewespecializeto sion for the average total current and for the noise one-dimensionalsystemsandinvestigatepumpingofelec- power of electrons pumped in a tight-binding wire. Al- tronsthroughthree different structures: a single barrier, ternatively, one can combine Floquet theory with non- aseriesofbarriersandaquantumwell. InSectionIVwe equilibrium Green’s function techniques.8,9 Generally summarize the main results and discuss future projects. 2 II. TIME-DEPENDENT CURRENT provided that a) t goes to infinity and b) the retarded Green’s function projected on the central region, GR, We consider an open quantum system C (central re- [or the advanced one, GA] vanishes when the separation gion) connected to two macroscopically large reservoirs betweenitstimeargumentsgoestoinfinity. Intheabove L and R (left and rightelectrodes). We are interestedin equation Σ = ΣL+ΣR is the embedding self-energy in describing the electron dynamics when region C is dis- the long time limit and the symbol Tr denotes a trace turbed by arbitrary time-dependent electric fields. As- over a complete set of states of the central region. We suming that the reservoirs are not directly connected, also have used the short-hand notation f g (t1;t2) the one-particle Hamiltonian of the entire system reads −∞∞dt¯f(t1;t¯)g(t¯;t2)fortheconvolutiono{ft·wo}function≡s fR and g. HLL HLC 0 For an applied bias Uα in reservoir α = L, R, which H(t)=HCL HCC(t) HCR . (1) is constant in time, the embedding self-energy depends  0 HRC HRR  only on the difference between its time arguments. Let TheHamiltonianH ,α=L,R,aswellastheHamilto- i αα ΣR/A(ω)=Λ (ω) Γ (ω) (4) nianofthe centralregionHCC,areobtainedbyproject- α α ∓ 2 α ingthe fullHamiltonianH ontothe subspaceofthe cor- be the Fourier transform of the retarded/advanced self- respondingregion. Howto choosethe one-particlestates energy. The imaginary part Γ is the contribution of in regions L, R or C depends on the specific problem α regionα tothe localspectraldensity. TheFouriertrans- at hand. We can use, e.g., a real-space grid for ab-initio form of the lesser self-energy is then given by calculations, or a tight-binding representation for model calculations,orevendifferentbasisfunctionsfordifferent Σ<(ω)=if (ω)Γ (ω), (5) regions (for instance, eigenfunctions of the reservoirs for α α α L and R, and localized states for C).10 The off-diagonal where f (ω)=f(ω U ) is the Fermi distribution func- α α − partsinEq. (1)accountforthecontactsandaregivenin tion. terms of matrix elements of H between states of C and Letus specialize to periodic time dependent perturba- states of L and R. tions in region C: H (t) = H (t+T ). According CC CC 0 In many applications of physical interest the driving to Floquet theory, we assume that the Green’s function fieldisperiodicintime. Inthiscaseitispossibletowork in Eq. (3) can be expanded as follows out an analytic expression for the dc component of the total current, Idc, provided memory effects and initial- GR/A(t;t′)= dωGR/A(ω)e−iω(t−t′)+imω0t′, (6) state dependence are washed out in the long time limit. XZ 2π m m Below we combine the Floquet formalism with nonequi- where ω = 2π/T is the frequency of the driving field. librium Green’s functions and generalize the formula for 0 0 I by Camalet et al.7 to arbitrary contacts. We also We wish to emphasize that the above expansion is justi- dc fied only if all obervable quantities (calculable from G) discuss the limitations of Floquet theory and propose an oscillate in time with the same frequency as the external alternative approachbased on the real time propagation field. As pointed out by Hone et al.,16 this is a question- of the initially occupied states of the system. able assumption. Inserting Eq. (6) into Eq. (3) and extracting the dc A. Long time limit: Floquet theory and Keldysh component we obtain formalism 1 t+T0 Q lim dt¯Q (t¯) Mostapproachestodrivennanoscalesystemsarebased α,dc ≡ t→∞T0 Zt α onafictitious partitioningfirstintroducedby Caroliand = dωG (ω)Σ<(ω)+ dω coworkers.11 The initial many-particle state is a Slater Z 2π 0 α Z 2π determinant of eigenstates of the isolated left and right G (ω)Σ<(ω)G† (ω)ΣA(ω mω ),(7) reservoirs with eigenenergy below some chemical poten- × X m m α − 0 m tial µ. A more physical initial state has been considered by Cini.12 It is a Slater determinant of eigenstates where we have defined of the contacted system L + C + R with eigenenergy G (ω) GR(ω mω )=[GA (ω)]†. (8) smaller than µ. Independently of the initial state, it has m ≡ m − 0 −m been proved13,14 that the number of electrons per unit The last equality in Eq. (8) follows directly from the time that leave the α = L, R reservoir is given by the identityGR(t;t′)=[GA(t′;t)]†. The dccomponentI formula15 α,dc of the time dependent total current I (t) is givenby the α righthandside ofEq. (2)with Q (t) replacedby Q . Iα(t)=2ReTr[Qα(t)], (2) α α,dc InAppendix Aweshowthatinthemonochromaticcase, Q (t)= GR Σ<+GR Σ< GA ΣA (t;t), (3) H (t)=H0 +U eiω0t+U e−iω0t, (9) α { · α · · · α} CC CC + − 3 the resulting expression for I can be cast in a E = m ω is much larger than any other energy α,dc max max 0 Landauer-like formula scale in the problem. Typically, m is in the range max from 10 to 100. The coefficients G (ω), m > 0, are dω ±m I = [f (ω)T (ω) f (ω)T (ω)], then calculatedfromG (ω) by simple matrixmul- L,dc Z 2π L m,L − R m,R ±(m−1) Xm tiplications according to Eq. (A20). Knowing the Gm’s (10) one can compute the inelastic transmission probabilities and IR,dc = IL,dc, as it should be due to charge con- from Eqs. (11,12), and hence the averagecurrent. − servation. The “inelastic” transmission coefficients T m,α may be interpreted as the probability for electrons to be Intheaboveproceduremostofthecomputationaltime transmitted from one reservoir to the other with the ab- isspentformatrixinversionsandmatrixmultiplications. sorption/emissionofm quantaofthe drivingfield. They We canroughlyextimatethe overalltime ofafullrunas can be written as T m N (τ +τ ), where N is the number run max ω i m ω ≃ × × of mesh points (generally not uniform) along the ω axis T (ω)=Tr Γ (ω)G† (ω)Γ (ω mω )G (ω) , m,L h L m R − 0 m i used to evaluate the integral in Eq. (10), and τi (τm) (11) is the time for a single matrix inversion(multiplication). In our case both τ and τ scale as N3. Depending on i m Tm,R(ω)=TrhΓR(ω)G†m(ω)ΓL(ω−mω0)Gm(ω)i. thesystemandontheexternaldrivingforces,theinelas- (12) tic transmission probabilities might exhibit quite sharp WeobservethatforzerodrivingtheFouriercoefficients peaks as function of energy. Therefore for an accurate G , and hence the transmission probabilities T , are computation of the energy integral in Eq. (10) a fine m m,α all zero except for m = 0, and Eq. (10) reduces to the energy grid is required, which means that Nω is large. well-knownLandauerformulaforsteady-statecurrents.17 In the numerical calculations of Section III Nω is in the Ocunrrtehnetcwohnetnraaryd,railvlitnhgefiTemld,αi’sspcroensternibt.utTehteoctohrereasvpeorangde- rtoan1g0e51.00to1000. WeconcludethatTrun/(τi+τm)∼103 ing G ’s can be computed recursively from the zero-th m ordercoefficientG ,andwedeferthereadertoAppendix (ii) Periodic potentials. Beyond the monochromatic 0 Aforapracticalimplementationscheme. Itisalsoworth case, the recursive scheme of Appendix A becomes com- emphasizingthatourformulafortheT ’scorrectlyre- putationally demanding. The inclusion of one, two, m,α duces to the one of Camalet et al.7 for a central region ...more harmonics in the expansion of the driving field describedbyatight-bindingwireofsites 1 ,..., N and [see Eq. (A4)] transforms the block tridiagonal system | i | i connected to the left reservoir through 1 and to the of equations for the G ’s into a block penta-diagonal, m | i right reservoir through N . In this case, the spectral hepta-diagonal, ...system of equations. For arbitrary density matrices Γ hav|e oinly one nonvanishing entry, periodic drivings a Floquet-based approach may not be α [Γ ] = δ δ γ and [Γ ] = δ δ γ , and feasible. L n,m n,1 m,1 L R n,m n,N m,N R the coefficients T can be rewritten as m,α (iii) Arbitrary time dependent potentials. Besides the Tm,L(ω)=γL(ω)γR(ω mω0) [Gm(ω)]N,1 2, (13) wide class of periodic drivings, it is of interest to inves- − | | tigate the response ofa nanodeviceto non-periodic driv- T (ω)=γ (ω)γ (ω mω ) [G (ω)] 2. (14) ings as well.19 In such cases the Floquet formalism does m,R R L 0 m 1,N − | | notapplyandafulltime-dependentapproachisrequired. Equation(10)demonstrateshowtheinitialFloquetas- sumption of Eq. (6) allows for carrying the analytic cal- (iv) Transients. The Landauer formalism provides culationofthecurrent[Eq. (2)]muchfurtherandeventu- a very powerful technique to calculate non-equilibrium allydeliversa simple numericalschemefor the computa- quantities in steady-state regimes. Similarly, the Flo- tionoftheaveragecurrent. Despitetheenormoussuccess quetformalismallowstocalculatenon-equilibriumquan- in predicting ac dynamical properties of many different tities in “oscillating-state” regimes, i.e., when all tran- nanoscale conductors, Floquet theory might be inade- sient effects are died off. However, transient responses quate to face the future challenges of nanotechnology.18 can be expected to become of some relevance in the fu- Below we discuss some limitations of Floquet-based ap- ture. Molecular devices will eventually be integrated in proaches. nanoscale circuits and respond to ultrafast external sig- (i) Numerical performance. For later comparisonwith nals. Transienteffectsinsuchoperativeregimesmaynot ourproposedrealtimeapproach,webrieflyreportonthe be irrelevant, as it has been recently recognized by sev- numericalperformanceofFloquetalgorithms,likethere- eral authors.13,20,21,22,23,24,25 In Section III we provide cursive scheme proposed in Appendix A. Let N be the explicit evidence of long-lived superimposed oscillations number of basis functions in region C. For a given fre- inthetime-dependentcurrentprofile. Thefrequenciesof quency ω the calculationof G (ω) requires the inversion these oscillations are not commensurable with the driv- 0 ofm complexmatricesofdimensionN N. Thenum- ing frequency,andhaveto be ascribedto the presenceof max ber m should be chosen such that the×cut-off energy “adiabatic” bound states.26,27 max 4 B. Real time propagation to second-order in δ. From Eq. (18) we can extract an equationforthetime-evolvedstateinregionC,similarly In this Section we propose an alternative approch to towhatwehavedoneforthederivationofEq. (16). The driven nanoscale transport. The main idea is to cal- final result is culate the time-dependent total current from the time- 1 +iδH(m) ψ(m+1) = 1 iδH(m) ψ(m) dependent wavefunctions |ψs(t)i, where |ψs(0)i is the s- (cid:16) C eff (cid:17)| C i (cid:16) C − eff (cid:17)| C i th eigenstate of the system L+C +R before the time- + S(m) M(m) , (19) dependent perturbation is switched on. Our approach | i−| i doesnotrelyontheFloquetassumption,andisfreefrom with 1 the identity matrix in region C. Equation (19) C allthelimitationsdiscussedpreviously. Furthermore,the is the proper (unitary) time-discretization of Eq. (16). computationaltimeiscomparablewithFloquet-basedal- Moreover, Eq. (19) is ready to be implemented since it gorithms. contains only finite-size matrices and vectors (with the As the full Hamiltonian H(t) refers to an extended dimensionusedtodescribethecentralregionas,e.g.,the andnon-periodicsystem,wecannotsolvebrute force the number of lattice sites in a tight-binding representation Schrödinger equation or the number of grid points in a real-space grid representation). In the following we give full implementation d i ψ(t) =H(t)ψ(t) . (15) details of the various terms in Eq. (19). dt| i | i Forthesakeofsimplicity,weconsiderone-dimensional semi-infinite reservoirs described by tridiagonal matri- Fortunately,wedo notneedto calculatethe time depen- ces H , α = L, R, with diagonal entries h and off- dent wavefunction everywhere in the system in order to αα α diagonalentriesV ,seeFig. 1. Fortight-bindingmodels, calculate the total current. The knowledge of the wave- α the parameterh representsthe on-site energywhile the function in region C is enough for our purposes (see be- α parameter V represents the hopping integral between low). Denoting with ψ (t) the wavefunction projected α C | i nearest-neighbour sites. The Hamiltonian H is also onregionC andwith ψ (t) thewavefunctionprojected αα α | i suited to describe continuum models with a three-point on region α = L, R, it is straightforward to show that Eq. (15) implies the following equation for ψ (t) 28 discretization of the kinetic term. In this case, the pa- | C i rameter h = 1/∆x2+U and V = 1/(2∆x2), where α α α − d t ∆x is the gridspacing. We wouldlike to emphasize that idt|ψC(t)i = HCC(t)|ψC(t)i+Z ΣR(t;t′)|ψC(t′)i thealgorithmcaneasilybegeneralizedtoreservoirswith 0 anarbitrarysemi-infiniteperiodicityanditisnotlimited + HCαexp( iHααt)ψα(0) , (16) to one-dimensional systems.21 − | i α=XL,R Withoutlossofgenerality,weconsideracentralregion thatincludesfewsitesoftheleftandrightreservoirs,and where we denote by α the state where only the site of region | i C connected to the reservoir α = L, R is occupied (see ΣR(t;t′)= H exp( iH (t t′))H (17) Cα − αα − αC Fig. 1). X α=L,R Region C is the Fourier transform of the embedding self-energy in Eq. (4). ... hL hL hR hR ... Equation (16) is an exact equation for the time evo- V V V V L L R R lution of open systems, but is still not suited for a nu- |2;L> |1;L> |L> |R> |1;R> |2;R> merical implementation. The importance of charge con- servation in quantum transport makes the unitary prop- FIG. 1: The schematic sketch of the electrode-junction- erty a fundamental requirement. In this work we use electrode system with semiperiodic one-dimensional elec- a unitary algorithm which has been recently proposed trodes. to study electron transport in biased electrode-device- electrode systems.21 Below we illustrate the main ideas The memory state M(m) stems from the the second andspecializetheformulasofRef. 21toone-dimensional | i term on the r.h.s. of Eq. (16) and reads reservoirs. For a given initial state ψ(0) = ψ(0) we calculate M(0) =0, (20) the time-evolved state ψ(t| = i2mδ)| =iψ(m) by ap- | i m | i | i proximating Eq. (15) with the Crank-Nicholsonformula while for m 1 we have ≥ 1+iδH(m) ψ(m+1) = 1 iδH(m) ψ(m) , (18) m−1 (cid:16) (cid:17)| i (cid:16) − (cid:17)| i M(m) = δ2 α αψ(k+1) + αψ(k) | i X | i X hh | C i h | C ii α=L,R k=0 with H(m) = 1[H(t )+H(t )]. The above propa- gationschemei2sunitmar+y1(normcomnserving)andaccurate × hqα(m−k)+qα(m−k−1)i. (21) 5 Theq-coefficientscanbecomputedrecursivelyaccording Finally, the effective Hamiltonian is given by to H(m) =H(m) iδ q(0) α α. (31) (1+iδh )+ (1+iδh )2+(2δV )2 eff CC − X α | ih | q(0) = − α α α , (22) α=L,R α p2δ2 The above algorithm allows us to calculate the time evolution of any initial state whose wavefunction in the 1 iδh 2δ2q(0) reservoirs has the form in Eq. (26). This is the case of qα(1) = 1+−iδhα+−2δ2qα(0)qα(0), (23) both the contacting approach by Caroli and coworkers α α and the partition-free approach by Cini. In the former approach,theinitialone-particlestatesareeigenstatesof and for m 2 ≥ the isolated left and right reservoirs,meaning that qα(m) = qα(1)qqα(α(0m)−1) −δ21+iqδα(0h)αqα(+m−22δ)2qα(0) |ψα(0)i=2X∞ sin(pαj)|j;αi= |+pαi−i |−pαi, (32) m−1q(k)+2q(k−1)+q(k−2) j=1 δ2 α α α q(m−k), (24) − kX=1 1+iδhα+2δ2qα(0) α zfoerroαin=reLgi(oonrCα.=InRth),ezlaetrtoerfoarpαpr=oacRht(horecαom=pLu)t,atainodn oftheinitialone-particlestatesismoreinvolved. Herewe with the convention that q(m) =0 for negative m. α have used a recently proposed general scheme based on The source state S(m) stems from the last term on thediagonalizationoftheimaginarypartoftheretarded | i the r.h.s. of Eq. (16) and reads Green’s function.21 This scheme may also be used for arbitrary, semiperiodic electrodes. In the special case of (1 iδH )m S(m) = 2iδ H α− αα ψ (0) , spatially uniform one-dimensional reservoirs one can, of | i − X Cα(1α+iδHαα)m+1| α i course, always use the textbook procedure of matching α=L,R (25) the wavefunction at the interfaces. where 1α is the unit matrix in region α. The source Denoting with |ψs,C(t)i the evolution of the original state depends on the initial wavefunction in the reser- eigenstate ψs(0) in the central region, we can calculate | i voirs. As we are interested in propagating eigenstates of the time dependent occupation ρ(j,t) of a state j in | i H(0), ψ (0) has the following general expression region C according to α | i |ψα(0)i=A(α+)|+pαi+Aα(−)|−pαi, (26) ρ(j,t)=Xf(εs)|hj|ψs,C(t)i|2, (33) s with where ε is the eigenvalue of ψ (0) and f(ε) is the s s | i ∞ Fermi distribution function. Similarly, the total time- p = eipαj j;α , (27) dependentcurrentI (t)canbecalculatedfromthetime- α α | i | i Xj=1 derivative of the total number of particles in electrode α and reads and the state j;α where only the j-th site of reservoir | i α = L, R is occupied, see Fig. 1. For extendend states Iα(t) = 2 f(εs) − inregionαthe parameterpα isreal. Forboundstatesor Xs jX6=α fullyreflectedstatesinregionαtheparameterpαisimag- Im j ψs,C(t) ψs,C(t)α [HCC(t)]αj, (34) inary and the amplitude (A(+) or A(−)) of the growing × h | ih | i α α exponentialis zero. No matter if p is realor imaginary, where the sum is over all states j of region C except the α one can prove that state α . | i We wish to conclude this Section with a discussion on (1 iδH )m the performance of our method and a comparison with H α− αα p =ζ(m) α , (28) Cα(1 +iδH )m+1| αi α | i Floquet-based approaches. α αα (i) The computational time T scales linearly with run with the number of states N used to evaluate the sum in Eq. s (33) or Eq. (34), and quadratically with the number of m time steps N . In most cases transient effects disappear ζ(m) =eipαV γ(m)+iδ γ(m−k) q(k)+q(k−1) , (29) t α α α X α h α α i after few femtoseconds (few tens of atomic units). Using k=0 a time step of the order of 10−2 a.u. we can obtain a rathergoodestimateofI withN 103to104. Given and α,dc t ∼ a central regions with hundreds of states the real time γ(m) = (1−iδhα−2iδVαcospα)m . (30) algorithmcanbeofthesamespeedoforevenfasterthan α (1+iδh +2iδV cosp )m+1 the Floquet algorithm of Appendix A. α α α 6 (ii) The real-time algorithm can deal with arbitrary spatiallyrestrictedtotheexplicitlytreateddeviceregion (periodic and non-periodic) drivings, and the computa- which in our case also coincides with the static potential tional time is independent of the specific time depen- barrier. The barrier extends from x = 8 to x = +8 − dence of H (t). Moreover, the algorithm is easily a.u. and its height is 0.5 a.u., see inset (b) in Fig 2. The CC generalized21 to deal with spatially uniform bias poten- system is unbiased, i.e., U = U = 0, and the Fermi L R tialsinthe electrodeswitharbitrarydependenceontime energy of the initial (ground) state is ε = 0.3 a.u.. For F such as, e.g., for an ac bias. the numerical implementation we have chosen a lattice (iii) From the time-evolved states ψ (t) we have ac- spacing ∆x=0.08 a.u., and 200 k-points between 0 and s | i cesstothetotalcurrentI (t)atanytimet,andnotonly k = √2ε which amounts to the propagation of 400 α F F to the long-time limit of the dc component of I (t). In states. α particular, we can easily investigate transients and the full shape of Iα(t) for t . In practice, this limit is 0.02 12 < I ( t ) > ( 10 - 4 a.u.) (a) → ∞ achieved for a finite time T after which all transient max 8 left phenomena have died out. 0.01 middle niidntaeyip(ltiiepovanf)oldtisAeentnnancttltouieatddhleiiseninsrstgilhatoeeydclvaegfaulcrntoinnrtucoatnbngidoo-eentslhoaetflacsttotprehuaoerconoefmriatyneh.n2tted9ehrocItaondicmndtiteiesaoepcntdthrs,eoetdvvhpiiesadoyeesstdsxtiimtebtemhierle---. < I ( t ) > ( a.u.) 0 4050.1 100 1 t (a.5u0.0) 10 100right t( b(a).u.) Therefore, according to the Runge-Gross theorem,29,30 -0.01 the interacting time-dependent density can be repro- duced in a fictitious system of non-interacting electrons -0.02 x = -8 x = 0 x = 8 moving under the influence of an effective Kohn-Sham potential which is local in space and time. We observe thatthis isnotthe caseinthe contactingapproachsince FIG.2: Time-dependentaveragecurrentattheleftandright theswitchingofthecontactsmakestheexternalpotential interfaceandinthemiddleofregionC forpumpingthrougha non-local in space and hence the Runge-Gross theorem single squarebarrier bya travelling wave. The travelling po- does not apply. tential wave is restricted to the propagation window |x| < 8 (v) Finally, we would like to stress that the Hamilto- a.u. and has the form U(x,t) = U1sin(qx−ω0t) with am- nianH (t)entersinthealgorithmonlyviatheeffective plitude U1 = 0.35 a.u., wave number q = 1.6 a.u. and fre- CC quencyω0 =0.2a.u.. Inset(a)isamagnificationofthelong- Hamiltonian H of Eq. (31), and has no restrictions. eff time behavior. The straight line corresponds to the value Thus,besidesone-dimensionalstructures(likethosecon- IL,dc = 7.63 ·10−4 a.u. of the average current calculated sidered in Section III) one can consider other geometries usingtheFloquetalgorithm. Inset(b)displaysthestaticpo- as well, like those of planar molecules, quantum rings, tentialbarrier(solidline)andthesuperimposedright-moving nanotubes, jellium slabs, etc. travelling wave (dashed line). In Fig. 2 we plot the time-dependent average current III. NUMERICAL RESULTS calculated according to In this Section we illustrate the performance of the 1 t I(t) = θ(T t) dτI(τ) proposed scheme by presenting our results for one- h i 0− t Z 0 dimensional continuos systems described by the time- 1 t dependent Hamiltonian + θ(t T ) dτI(τ), (36) − 0 T Z 0 t−T0 2 H(x,t)= ∇ +U(x,t). (35) with the period of the travelling wave T = 2π/ω . For − 2 0 0 the time propagation we have chosen a time step 2δ = We discretize H on a equidistant grid and use a three- 0.02 a.u.. As expected I(t) converges to some steady h i point discretization for the kinetic term. Within this value I after a transient time of the order of 50 60 L,dc ÷ model we study various model systems highlighting dif- a.u.. We have calculated the average current in three ferent features in electron pumping. different points of region C and verified that the steady value does not depend on the position. The dc limit I can also be computed from the Floquet algorithm L,dc A. Archimedean screw of Appendix A. Using m =15 and N =150 energy max ω points between0andε , wefind I =7.63 10−4 a.u., F L,dc · As a first example of electron pumping we have cal- in very good agreement with the average current of the culated the time evolution of the density and total cur- time propagated system, see inset (a) of Fig. 2. rent for a single square barrier exposed to a travelling In Fig. 3 we plot the time-dependent density ρ(x,t) potential wave U(x,t) = U sin(qx ω t). The wave is in the device region as a function of both position x and 1 0 − 7 200 inversion. The present Section is not intended to give a me t (a.u.1)00 150 realistic description of some specific experiment. ti Floquet 50 left 0.1 middle ) 0.06 (a.u. right ).u.a( ytisned 0.04 < I ( t ) > 0 0.1 1 10 100 t (a.u.) 0.02 -0.1 t > 0 -8 -4 0 position x (a.u.) 4 8 0 t < 0 -0.2 x = -6 x = 0 x = 6 FIG. 3: Time-dependent density in region C as a function of position x and time t. The range for ρ(x,t) is between 0 and FIG.4: Time-dependentaveragecurrentattheleftandright 0.06a.u. forimprovedvisibilityofthepocketsofdensity. All interface and in themiddle of region C for pumpingthrough parameters are thesame as for Fig. 2. a device region which extends from x = −6 to x = 6 a.u.. A travelling wave U(x,t) = U1sin(qx−ω0t) with U1 = 0.5 a.u., ω0 = 0.8 a.u., and q = 0.6 a.u., is superimposed to the time t. The density exhibits local maxima in the poten- static potential U0(x)=UC(1+cos(kx)) with UC =0.5 a.u. and k = 10π/6 ∼ 5.2 a.u., and the all system is unbiased, tial minima and is transported in pockets by the wave. see inset. The straight line corresponds to the value of the From Fig. 3 it is also evident that the height of the averagecurrentasobtainedfromtheFloquetalgorithmwhich pockets is not uniform over the system, and reaches its yields IL,dc=3.26·10−2 a.u.. maximum around x = 0. We also notice that the particle current flows in the same direction as the driving In Fig. 4 we plot the time-dependent average current wave. The pumping mechanism in this example resem- [see Eq. (36)] in three different points of the device re- bles pumping of water with the Archimedean screw. gion. For the numerical propagation we have used a lattice spacing ∆x = 0.06 a.u., a time step 2δ = 0.02 a.u., and 400 k-points between 0 and k = √2ε . The sys- F F B. Pumping through a semiconductor tem responds to a right-moving travelling wave by gen- nanostructure erating a net current flowing to the left. Again we observe that the transient time is of the order of few tens The second example was motivated by a recent exp- ofatomicunits,andthatthesteadyvalueisindependent eriment on pumping through a carbon nanotube.3 The of the position. As in the previous example, we used the arrangement has been suggested by Talyanskii et al.31 Floquet algorithm of Appendix A for benchmarking our and is as follows. A semiconducting nanotube lying on real-time propagation algorithm. Due to the high Fermi a quartz substrate is placed between two metallic con- energy the calculation was carried out with mmax = 15 tacts. A transducer generates an acoustic wave on the and Nω = 400 energy points between 0 and εF. The re- surface of the piezoelectric crystal. The crystalresponds sult IL,dc = 3.26 10−2 a.u. is displayed in Fig. 4 with · by generating an electrostatic potential wave which acts a straight line and is in extremely good agreement with likeourtravellingwaveontheelectronsinthenanotube. the long-time limit of the average current obtained from The directionofthe pumping currentis foundto depend direct propagation in time. on the applied gate voltage. A pumping current flowing To understand how the electron fluid moves when the in the direction opposite to the propagation direction of direction of the current is opposite to that of the driv- the travelling wave has been interpreted in a stationary ing potential wave, we have studied the dynamical flow picture as a predominant hole tunneling over electron pattern of the density. We emphasize that such a study tunneling. To reproduce such an inversion in the cur- wouldhavebeenrathercomplicatedinFloquet-basedap- rentflowwehavemodelledthe nanotube withaperiodic proaches. Thelatterareoftenusedas“blackboxes”,and static corrugation U (x)=U (1+cos(kx)) in region C, one needs to resort to limiting cases, like, e.g., the adi- 0 C with U = 0.5 a.u. and k = 10π/6 5.2 a.u. (see inset abatic picture, the high frequency limit, the theory of C ∼ inFig4). ForatravellingwaveU(x,t)=U sin(qx ω t), linear response, etc., for a qualitative understanding. 1 0 − with U = 0.5 a.u., ω = 0.8 a.u., and q = 0.6 a.u., we InFig. 5wedisplayacontourplotoftheexcessdensity 1 0 havefoundthattheminimum currentoccursatε =3.0 ∆ρ(x,t) = ρ(x,t) ρ(x,0) in an extended region which F − a.u.. All parameters in this example have been chosen includes the device regionandaportionofthe left reser- to better illustrate and discuss the effect of the current voir. In the device region we clearly see pockets that 8 1 R T L / T L(ε) T (ε) 0.5 R ε (a.u.) 0 0.03 u.) a. I ( L,dc 0 -0.03 ε (a.u.) F 0 1 2 3 4 5 6 FIG. 6: dc component of pump current IL,dc and left/right transmission probabilities T as function of the Fermi en- L/R ergy. ThecurveshavebeenobtainedusingtheFlocquetalgo- FIG. 5: Contour plot of the excess density ∆ρ(x,t) in the rithm of Appendix A with mmax =15 and Nω =800 energy device region (x between -6 and 6 a.u.) and in a portion of points between 0 and εF=6 a.u.. theleftreservoir(xbetween-30and-6a.u.). Duetothelarge oscillationsoftheexcessdensityinthedeviceregion,∆ρ(x,t) has been scaled down by a factor of 10 for |x| < 6 a.u.. We abilities T = T , α = L,R, for left- and right- also draw straight lines to show a pocket trajectory and the α m m,α going electronsP[see Eqs. (11, 12)]. As one can see from trajectory of a superimposed oscillation. All parameters are Fig. 6, both T and T remain quite small for Fermi thesame as for Fig. 4. L R energiesbelow∆ 0.54a.u.,whichroughlycorresponds ∼ to the bottom of the lowest band of the periodic structure of the device. In this energy window transport is are dragged by the travelling wave and are moving to dominated by tunneling and the pumping current fol- the right. However,every pocket with a slightly positive lows the travellingwave(T >T ) similar to the case of L R ∆ρ is followed by a pocket with a noticeably negative the Archimedean screw, see Section IIIA. For ε > ∆ F ∆ρ, and the net excess density is negative. On the other we enter the region of resonant transport (the energy of hand,intheleftreservoironlypocketswithapositive∆ρ the lowest band) and T , T sharply increase. We ob- L R exist and they move to the left. We conclude that the servethat forε <ω =0.8 a.u. both T andT have a F 0 L R driving produces right-moving “bubbles” in the device structuresimilartothetotaltransmissionfunctionofthe regionandthattoeachbubblecorrespondsamoredense staticcase. Forω <ε <ω +∆,however,T decreases 0 F 0 L region of fluid moving to the left. One can estimate the significantly while T remains fairly constant around 1. R speedvpoc ofthe travellingpocketsfromthe slope ofthe We interpret this in the following way. The probability patternsatconstantdensityandfindvpoc ω0/q 1.33 of the right-going electrons of emitting a photon of fre- ∼ ∼ a.u., as expected. We also notice superimposed density quencyω (andthereforereducingtheir energy)is larger 0 oscillations on each pocket. These oscillations have the thanfortheleft-goingelectrons. Loosingthisenergy,the samespatialperiodofthestaticcorrugationinthedevice transmission T resembles the static transmission func- L region, and move in the same direction of the pockets at tion for energy ε ω which has a much lower value. F 0 a constant speed vosc ω0/k 0.15 a.u.. The asymmetry bet−ween left- and right-going states can ∼ ∼ InFig. 6 we illustrate howthe pumped currentinthis easily be understood by realizing that the pump wave model depends on the Fermi level. For Fermi energies introduces a preferential direction in the problem. As comparable to the amplitude of the corrugated poten- further evidence to support this interpretation we note tial in the device region the pumping current is always that for ε = ω +∆ the transmission function T in- F 0 L positive, i.e., follows the propagation of the perturbed creases rapidly as for ε = ∆. This can be viewed as F wave. However, there are striking effects that are more a replica of the static transmission function shifted by or less independent of the strength of the perurbation: one quantum of energy ω . Throughout the energy win- 0 the pumping current reaches a maximum positive value dow of the lowest band, T remains lower than T . As L R at ε ω = 0.8 a.u., then decreases with increasing a consequence the pumping current decreases monotoni- F 0 ∼ Fermi energy (with the turning point to negative val- cally. This behaviour is reversed when the Fermi energy ues just below ε = 2 a.u.) and reaches a minimum hits the top of the lowest band, around 3.4 a.u.. In the F (negative) value above ε = 3 a.u.. To rationalise this gap (of about 2U = 1 a.u.) both T and T drop and F C L R behaviorwehavecalculatedthe totaltransmissionprob- transportisdominatedbytunnelingagain. Inthisregion 9 T >T and the pumping current increases. 2δ =0.05a.u.,andcalculatethecurrentI(t)atthecenter L R The present model gives positive and negative pump- ofthequantumwell. AsintheexamplesofSectionsIIIA ingcurrentasafunctionoftheFermienergyandprovides and IIIB, one observes a first transient behavior which a simple physical interpretation of the effect of current lasts for few tens of atomic units. However, after this inversion. Our picture, however, is somewhat different first normal transient a second transient regime sets in. from the one given by Leek et al..3 Indeed, in their ex- In Fig. 7 we plot the modulus of the discrete Fourier planationthesignofthepumpingcurrentisindependent transform of the current of the frequency ω of the travelling wave. On the other 0 hand, in our case, if the frequency exceeds the width of np+N0 2πk the lowest band, the right-going electrons cannot emit a I˜(ωk)=2δ I(2nδ)e−iωknδ, ωk = (37) N δ X 0 photonandcurrentinversionisnotguaranteedanymore. n=np for n =(4+2p) 103, p=0,1,2,3,4,and N =16 103 p 0 · · C. Transients effects (correspondingtothetime intervalst (t ,t +T )with p p 0 ∈ t = (2+p) 100 a.u. and T = 800 a.u.). Besides the p 0 × As alastexamplewestudy electronpumpinginquan- zero-frequency peak (not shown) due to the non vanish- tum wells. We will show the presence of long-lived su- ingdccurrent,thestructureofI˜(ω)hasfivemorepeaks. perimposed oscillations whose frequency is generally not Belowwe discussthe physicaloriginofthese extrapeaks commensurable with the driving frequency. The quan- and show that they are related to different kinds of in- tumwellismodelledwithastaticpotentialU (x)= 1.4 ternal transitions. 0 − a.u. for x < 1.2 a.u. and zero otherwise. Initially Wefirstobservethatthebiasedsystemhastwobound the system| |is in the ground state with Fermi energy states with energy ε∞ = 1.032 a.u. and ε∞ = 0.133 b,1 − b,2 − εF = 0.1 a.u.. The unperturbed Hamiltonian has two a.u. (slightly different from the bound-state energies of bound-state eigensolutions with energy ε0 = 1.035 theunbiasedsystem). Thefirstandthelasttwopeaksoc- b,1 − a.u. and ε0 = 0.156 a.u.. The ground-state Slater cur at the same frequency of the bound-continuum tran- determinantb,2conta−insallextendedstateswithenergybe- sitions ε∞b,i → εF, and ε∞b,i → εF + UR, with i = 1,2. tween 0 and ε and two localized states with negative These sharp structures (mathematically stemming from F energy. At positive times a constant bias U = 0.1 thediscontinuityofthezero-temperatureFermidistribu- R a.u. is applied on the right lead and a travelling wave tion function) give rise to long-lived oscillations of the U(x,t)=U sin(qx ω t),withq =0.5a.u. andω =0.5 total current and density. Such an oscillatory transient 1 0 0 − a.u., is switched on in the quantum well. In the nu- regimediesoffslowlyas1/t. Thepower-lowbehaviorcan merical simulations we set the propagation window be- alsobe seen in the inset ofFig. 7, where a magnification tween x = 1.2 and x = 1.2 a.u. (which coincides with oftheregionwithtransitionsfromtheweaklyboundelec- − the static potential well) and choose a lattice spacing tron to the two continua is displayed. Denoting with R p ∆x = 0.024 a.u.. The occupied part of the continuum the product between the height of the second peak and spectrum is discretized with 100 k-points between 0 and the propagationtime t +T we have found R =26.305 p 0 2 k =√2ε . a.u., R = 26.307 a.u., and R = 26.328 a.u., which is F F 3 4 infairlygoodagreementwiththe expected1/tbehavior. )u.0.14 0.03 ε8b , 2 −ε8b,1 Tinhgertefoarned,tahpephriogahcthoefstzheeropeianksthdeeclrimeaistestwith inc.reOasn- (a.0.12 0.02 the cpontrary, the sharp peak at ω = ε∞ pε→∞ (∞bound- ω) p = 0 b,1− b,2 ~ (I 0.1 0.01 pp == 12 bTohuenodsctirlalantsiiotniosno)frtehmeabionusnudn-cbhoaunngdedtrawnitshitiionncrdeoasninotgdtpie. 0.08 p = 3 p = 4 off, in agreement with the findings of Refs. 26,27. We 0.06 0.25 0.3 0.35 0.4 emphasize that these latter oscillations are an intrinsic 0.04 ε −ε8 εF +UR −ε8b,2 property of the biased system and have nothing to do F b,2 with external drivings. 0.02 εF −ε8b,1 εF +UR −ε8b,1 Having discussed the behavior of the system which is 0.2 0.4 0.6 0.8 1 1.2 ω (a.u.) biased but not driven, we now study transient regimes in the biased and driven system, i.e., U = 0. Using the 1 6 samenumericalparametersasinthepreviousexamplewe FIG.7: ModulusofthediscreteFouriertransform ofthecur- evolvethe (non-interacting)many-body state fromt=0 rent for zero driving and U = 0.1 a.u.. The inset shows R a magnification of the region with bound-continuum transi- tot=1200a.u. withatimestep2δ =0.05a.u.. InFig. 8 tions. Different curvescorrespond todifferent time intervals. weplot the discrete Fouriertransformofthe currentcal- culated in the middle of the quantum well for different Letusfirstconsiderthebiasedsystemwithnodriving, amplitudes of the travelling wave U1 = 0.00, 0.01, 0.03 i.e., U = 0. We propagate the (non-interacting) many- a.u.. The time interval used to evaluate I˜(ω) is from 1 body state from t=0 to t=1400 a.u. using a time step t = 200 a.u. to t = 1200 a.u.. As expected, I˜(ω) has a 10 000...000444 (a) (b) driven quantum well has a very rich transient structure. u.)111...222 000...000333 000000......111111242424 Thisisduetothepresenceofboundstateswhichcansub- a. 000...111 ω) ( 111 000000......000000121212 000000......000000686868 stantially delay the development of the Floquet regime. ~ (I000...888 000...222555 000...333 000...333555 000...444 000000......000000242424 000...888555 000...999 000...999555 111 U = 0. IV. CONCLUSIONS AND OUTLOOKS 000...666 1 U = 0.01 1 first 000...444 U1 = 0.03 harmonic Time-dependentgate voltagescanbe usedtogenerate ε8b , 2 −ε8b,1 a net current between unbiased electrodes in nanoscale 000...222 εF −ε8b,2 εF +UR −ε8b,2 hsaermcoonndic junctions. Most works focus on periodic drivings for which Floquet-based approaches provide a powerful ma- 000...222 000...444 000...666 000...888ω (a.u.)111 chinery to investigate the long-time behavior of the system. Combining Floquet theory with nonequilibrium Green’s functions techniques we obtained a general for- FIG. 8: Modulus of the discrete Fourier transform of the mulafortheaveragecurrentofmonochromaticallydriven current for the biased quantum well (U = 0.1 a.u.) per- R turbedby thetravelling wave U(x,t)=U1sin(qx−ωt),with systems in terms of inelastic transmission probabilities. q = 0.5 a.u. and ω = 0.5 a.u.. Three different amplitudes The case of polychromatic drivings, which has received U1 = 0.00, 0.01, 0.03 a.u. are considered. Inset (a) displays scarce attention so far, is analytically more complicated a magnification of the region with bound-continuum transi- and computationally rather costly. tions. Inset (b) shows a magnification of the region with the In this work we proposed an alternative approach bound-boundtransition and the second-harmonic peak. which can deal with monochromatic, polychromatic and nonperiodic drivings. The computational cost is inde- pendentofthe particulartime dependence ofthe driving wellpronouncedpeakatthedrivingfrequency(firsthar- potential. As an extra bonus we can investigate how the monic). Increasing the amplitude of the driving field the transient behavior depends on the initial state and on height of the first-harmonic peak increases and higher- the details of the switching process. The basic idea is to order harmonic peaks become visible (breakdown of lin- calculate the time-dependent density and current from ear response theory). This is clearly shown in inset (b) the time-evolved (non-interacting) many-particle state. where the second-harmonic peak is visible for U1 =0.03 This amounts to solving a single-particle Schrödinger a.u. butnotforU1 =0.01a.u.. ThestructureofI˜(ω)has equationforeachoccupiedeigenstateoftheunperturbed alsootherpeaksatfrequencieswhicharenotcommensu- system. We have given full implementation details of rable with the driving frequency. Such peaks are due to the time-propagationalgorithmand discussedits perfor- the presence of bound states in the biased-only system. mance. The generalization to two- or three-dimensional In inset (a) we show a magnification of the region with reservoirs can be worked out following the general lines bound-continuumtransitions. Thedrivingfieldbroadens of Ref. 21 and its implementation is in progress. thepeak-structure,thusspeedingupthepower-lawtran- We illustrated our scheme in one-dimensional struc- sientregime. Theshapeofthebound-boundtransitionis tures. First we studied pumping through a single bar- displayed in inset (b). The height of the peak decreases rier, and showed that the electrons are dragged by the with increasing amplitudes and the transition changes travelling wave and move in pockets. Second we studied from an infinitely long-lived excitation to an excitation pumping in semiconducting structures, and investigated with a finite life time. Let m be the smallest integer s the phenomenon of current inversion. In both examples forwhichm ω >max(ε∞ , ε∞ ); forsmallamplitudes s 0 | b,1| | b,2| the Floquet algorithm of Appendix A is used for bench- the life time is proportional to (1/U2)ms according to marking the long-time limit of the real-time simulations 1 the following reasoning. The retarded Green’s function and we have found an excellent agreement between the in region C can be written in terms of the embedding twoapproaches. Finally,weconsideredpumpingthrough self-energy of Eq. (4) and the Floquet self-energyΣR of a quantum well connected to biased reservoirs. The aim ac Eq. (A18). The Floquet self-energy generates replica of of this latter example is to show the existence of a long- thecontinuousspectrumwhichareshiftedbymultiplein- lived transient regime in rather common physical sys- tegersofω0 andcontributes tothe imaginarypartofthe tems. The transient oscillations are explained in terms Green’s function, GR. The leading-order contribution of bound-bound transitions and bound-continuum tran- of the m-th replica to ImGR scales like (U2)m. There- sitions. Theseoscillationsusuallyhavefrequencieswhich 1 fore,bound-statesimplepolesofGRgetembeddedinthe are not commensurable with the driving frequency and continuumspectrumofsomeofthereplicaandaquirean arethereforenotdescribedbytheinitialFloquetassump- imaginary partproportionalto (U2)m, with m the order tion. 1 ofthe replica. The leading-ordercontributiontothe life- The present work opens the path towards systematic time of the bound-bound excitation is then proportional studies of nanoscale devices as it is not restricted to to (1/U2)ms. linear response theory and can cope with general time- 1 In conclusion, we have shown that the biased and dependentaswellasspatialperturbations. Ourapproach