Electron Pumping under Non-Markovian Dissipation: The Role of the Self-Consistent Field

Focusing on electron transport through a periodically driven resonant tunneling diode, we study the generation of a non-vanishing dc-current by applying symmetry breaking external ac ﬁ elds with phase di ﬀ erence φ in a statically unbiased system. The e ﬀ ect of an environment is investigated using the system-bath Hamiltonian represented by the electron system coupled to harmonic oscillator modes with a Drude – Lorentz spectral density. To carry out simulations, we use the hierarchal equations of motion approach in the Wigner representation including a self-consistently constructed electric ﬁ eld that is determined from the electron distribution using the Poisson equation. We show that the maximal pumping current at a phase di ﬀ erence near φ = π = 2 is strongly in ﬂ uenced by the system-bath coupling strength. The e ﬀ ect of dissipation is diminished if the self-consistent part of the potential is ignored.


Introduction
Electron pumping, i.e., generating non-vanishing average dc electronic currents in statically unbiased systems by external ac-fields has been realized in a wide range of experiments, e.g., with quantum dots, [1][2][3] nanotubes, 4) semiconductor heterostructures, 5,6) and a Josephson junction array.7) Theoretical descriptions of the pumping effect using the scattering matrix are given in Refs. 8-10 wi10][11] In Refs.9 and 10 a double barrier resonant tunneling system is studied, and, in Ref. 11, the pumping effect for the free particle under a bi-harmonic driving (harmonic mixing) is discussed.
An early example of a similar strategy is the heuristic approach to investigating the effect of time-periodic driving on the current voltage characteristics of superconductorinsulator-superconductor junctions, given by Tien and Gordon. 12)7][18][19] The authors of Refs.13-17 consider tight-binding Hamiltonians, while those of Ref. 18 study an interacting two-level system.It has been stressed that coherent quantum pumping occurs because of the interference of energetically different transport pathways. 10) Itis thus intriguing to investigate the influence of a dissipative environment on this effect.20,21) For a two-site molecular wire, in Ref. 20, it is shown that, in the harmonic mixing case, phonon damping significantly increases the increase of the pumping current for certain parameters while decreases it for other case (shift from a sine-like curve to a cosine-like curve as a function of the parameter).The authors of this groundbreaking study used the Floquet picture in addition to an approximate Hartree-Fock decoupling scheme, in order to treat the coupling between the electron system and the vibrational modes.The current increase by coupling to phonon modes has been corroborated recently in a study of dissipative transport through a Cooper pair sluice using a non-Markovian equation of motion approach for weak system bath coupling strength.21) In the present study, however, to incorporate the coupling to the phonon environment, we use a numerically rigorous time-dependent propagation approach (no weak coupling assumption), the reduced hierarchal equations of motion (HEOM) method, reviewed in Refs.22-24.The HEOM method is a non-perturbative approach that can be converged asymptotically to the desired accuracy even under strong time-dependent perturbations at finite temperatures by increasing the number of hierarchal elements.24) The quantum suppression of ratchet rectification was studied in Wigner phase space using the HEOM formalism.25) In the electron transport context with open boundaries, the Wigner phase space representation of HEOM is appropriate, because this can handle inflow and outflow boundary conditions.In the pure quantum case (without coupling to an environment), the merits of the phase space description were discussed by Frensley.26,27) A previous application of the HEOM method to electron transport in a resonant tunneling diode, including dissipative system bath coupling was given by Sakurai and Tanimura.28,29) The paper is organized as follows: In the second section, the model Hamilton for resonant tunneling under timedependent gate fields is introduced. An nalytical formula for the adiabatic current is reviewed.Moreover, numerical results for the pumped current as a function of the phase difference between the applied driving fields in the adiabatic limit are given.Then, in Sect. 3, we riefly review the Wigner function approach to the electron transport through a system with open boundaries and its extension to the dissipative case using the HEOM methodology to allow for the description of non-Markovian system dynamics.At the end of Sect.3, numerical details are discussed.In Sect.4, numerical results in the time-domain for the average pumped current as a function of the phase difference between the applied gate fields in the case of a dissipative environment are presented.Conclusions and an outlook are given in Sect. 5.

Model Hamiltonian and Adiabatic Electron Pumping without Dissipation
We start by introducing the model Hamiltonian for resonant tunneling and reviewing the basic mechanism of tunneling transport through double barrier structures driven by an external field that breaks the time reversal symmetry, leading to the phenomenon of electron pumping.

Double barrier resonant tunneling
We will extend the investigations of resonant tunneling of electrons through a double barrier heterostructure U static , used in previous publications to study the current-voltage characteristics of a resonant tunneling diode without 26,27) and with 28,29) coupling to a Caldeira-Leggett bath.The most intriguing finding in previous work was the numerical reproduction of the existence of a region of negative differential resistance in the IðVÞ curve of the current as a function of bias voltage, [26][27][28][29][30][31] as well as the observation of hysteresis, [28][29][30] self-excited current oscillations, 28,29,31) and tristability 32) in this critical region.
In contrast to those previous studies, here, we consider the case of zero external bias, but instead we apply a timedependent electric field to the statically unbiased potential employed there.As in previous work, across the device depicted in Fig. 1, the effective electron mass is assumed to be constant at a value of m ¼ 0:067 m 0 , with the bare electron mass m 0 .The barriers as well as the spacer layers have a width of L br ¼ 2:825 nm and the barrier height is 0.27 eV, whereas the quantum well has a width of L qw ¼ 4:520 nm and the contacts have length L c ¼ 16:95 nm, leading to a total device length of L d ¼ 49:72 nm.The full potential is expressed as Here U static ðqÞ is the static double barrier potential depicted in panel (b) of Fig. 1 and U self ðq; tÞ ¼ ÀeÈðq; tÞ is the selfconsistent addition to the potential that has first been discussed in the present context by Kluksdahl et al. 30) It is obtained by solving the Poisson equation with the inhomogeneity evaluated from an equation of motion for the electron distribution n À ðq; tÞ as À @ 2 @q 2 ½Èðq; tÞ ¼ e½n þ ðqÞ À n À ðq; tÞ; ð2Þ with the dielectric constant ¼ 12:85 and an electron donor doping density in the contact regions [see panel (a) of Fig. 1] of n þ ¼ 2 Â 10 18 cm −3 .Furthermore, Vðq; tÞ is the timedependent external potential.It consists of two (different) time-dependent sinusoidal gate fields expressed as where V r ðtÞ ¼ U g sinð!tÞ; ð5Þ and L ¼ L br þ L qw .We stress that a (nontrivial) phase shift φ (≠ 0; 2) is introduced between the sinusoidal oscillations of the time-dependent barriers in order to break the timereversal symmetry. 10,33)Here the barriers oscillate in height.
In the Kramers-Henneberger frame, the barriers would oscillate laterally. 34)We note that, while a quantum ratchet system with biharmonic forces has been studied for periodic potential systems, [13][14][15][16][17][18][19]25) the present study focuses on a timedependent double barrier potential with inflow and outflow boundary conditions. Morever, here we include the effects of the self-consistent field on the electron dynamics to investigate a realistic nano-device situation.

Electron pumping
Before treating the above mentioned model in its full complexity and including additional coupling to the environment, first we review results for electron pumping in a simple resonant tunneling model.Simplifying the double barrier model given above is realized by considering the Hamiltonian with two δ-function barriers, separated by L, located symmetrically with respect to the origin.We note that the dimensionality of the pumping amplitude V j is energy times length, due to the delta function nature of the potential, and it should not be confused with the quantities V l;r introduced in Eqs. ( 4) and ( 5).The transmission in the time-independent case, V j ðtÞ ¼ V 0;j ¼ 0:27 eV Â 2:825 nm [¼ 0:54 a.u.(atomic units)], with barrier parameters chosen to match the value of barrier height times barrier width of the model in Ref. 29, consists of a series of resonant tunneling peaks, 35) and the first peak is displayed in Fig. 2. For reasons of comparison, we also depict the results for the corresponding finite barrier case 36) in this picture.
To obtain an impression of what to expect for our timedomain results to be presented below, we review some results for the adiabatic pumping current, i.e., the current for very slow external driving in the pure quantum case without dissipation.As shown in Ref. 10, the adiabatic current is given by Brouwer's formula 8) This result is correct (in the adiabatic limit) for arbitrary ratio of temperature to frequency but only gives a good approximation to the true current at finite frequencies for moderate driving strength. 10)The scattering matrix entering the expression for the current in the case of two oscillating delta function barriers with Hamiltonian ( 6) is given by where the abbreviations k ¼ ffiffiffiffiffiffiffiffiffi ffi 2mE p =ħ; p j ðtÞ ¼ V j ðtÞm=ħ 2 ; ðtÞ ¼ ½1 À ÁðtÞe ikL ; Þ=k have been used and the barriers, separated by L, are oscillating according to As a function of the phase difference ' 2 ¼ ' (setting ' 1 ¼ 0) and for different external pumping amplitude, V 1;1 ¼ V 1;2 , we plot the average current I ad,1 (flowing to the left) in Fig. 3.There we used ¼ 0:054 eV for the chemical potential, which is motivated by a previous resonant tunneling diode study. 37)Identical results are found for I ad,2 (flowing to the right) if we set For a frequency !¼ 4:13 Â 10 12 rad=s, corresponding to an energy smaller than the resonance width displayed in Fig. 2, and for relatively small driving strength, the adiabatic formula should be a good estimate for the true current. 10)For small driving amplitude the oscillation of the pumped current as a function of phase difference is almost sinusoidal with an extremum at ' ¼ =2, whereas for higher amplitude nonharmonic distortions of the current versus phase difference are observed and the extremal current is shifted towards smaller values of φ.We note in passing, that recently an investigation of three δ-function type barriers has been given, 38) which allows for studies of the pumped current as a function of two phase differences.

Wigner Phase Space Formulation of Dissipative Transport
In the following, we introduce a time-dependent viewpoint that allows us to describe the pumping effect not only under arbitrary driving frequency but also in presence of coupling to an environment.To this end, we first review the Wigner phase space approach briefly, and then we present the working formulae for the treatment of the non-Markovian system dynamics that arise due to the presence of a bath.

Wigner function description in the pure quantum case
The formulation of the open boundary conditions needed to describe the inflow of electrons from the boundaries and the outflow from the system is most straightforwardly achieved in Wigner phase space, as reviewed in the present context by Frensley. 27)rom the elements of the density matrix ðq; q 0 ; tÞ of a quantum system the corresponding Wigner function is calculated according to Wðp; q; tÞ ¼ 1 2ħ Z dr ðq þ r=2; q À r=2; tÞ exp½Àipr=ħ; ð11Þ with the usual factor of 1=2ħ included.The electron density in position space becomes n À ðq; tÞ ðq; q; tÞ ¼ Z dp Wðp; q; tÞ: The quantum Liouville equation for the Wigner function reads where the potential kernel in the quantum mechanical Liouvillian is given by We note that the nonlocal potential term in Eq. ( 13) can be written in two alternative forms, given by Groenewold 39) and by Moyal. 40)In the classical limit (ħ !0) this term becomes local and reduces to (using integration by parts) U W ðp À p 0 ; q; tÞWðp 0 ; q; tÞ !ħ!0 @Uðq; tÞ @q @W @p : The drift term [first term on the RHS of ( 13)] is of the same structure classically.The corresponding classical equation for the distribution function is the limit for vanishing damping strength of the Klein-Kramers, Fokker-Planck equation that has been discussed 41) and implemented 42) elsewhere.
The inflow conditions appropriate for the present situation are graphically depicted in Fig. 4, which amounts to setting Wðp; with the distribution function of the left and right reservoir, respectively.The particles leaving the device depend only on the state of the device, which is calculated by solving the Liouville equation and is not(!) fixed by the reservoirs.We note that there are conceptual problems of fixing the inflow boundary conditions in the quantum case without dissipation. 43)In the following we focus on the dissipative case.

Quantum hierarchal Fokker-Planck equations
The effect of localized surface vibrational modes on resonant tunneling has been investigated in several model inelastic tunneling studies. 44,45)The number of vibrational degrees of freedom that can be taken into account explicitly is rather low, however.In the following, we will allow for a continuous spectral density of the oscillators by moving to a reduced description of the dynamics.
As the starting point to describe the influence of environmental degrees of freedom on the electron pumping effect, we consider the Caldeira-Leggett Hamiltonian 46) for the electron with coordinate q coupled to bosonic degrees of freedom with m j , pj , xj , and !j being the mass, momentum, position and frequency of the jth phonon oscillator mode.Subsequently, we will restrict investigation to a bilinear coupling by setting Vð qÞ ¼ q.
The heat bath is characterized by its inverse temperature ¼ 1=kT and its spectral density which in the following will be assumed to be continuous and of Drude form 47) Jð!Þ ¼ ħm For factorized initial conditions and after tracing out the bath degrees of freedom in the path integral formalism, 48) the reduced quantum dynamics of the electron includes memory effects.The kernels of the time integrals that appear in the Feynman-Vernon influence functional are proportional to the canonical and symmetrized correlation function, respectively, of the collective bath coordinate For the Drude spectral density (20) they are given by 22,49) ÉðtÞ ¼ me Àjtj ; ð24Þ with the Matsubara frequencies k ¼ 2k=ħ 50) and Using the Wigner distribution and the quantum Liouvillian, the reduced equations of motion for the electron can be expressed in the form of quantum hierarchal Fokker-Planck (QHFP) equations in real time as @ @t W ðnÞ j 1 ;...;j K ðp; q; tÞ " # W ðnÞ j 1 ;...;j K ðp; q; tÞ þ È W ðnþ1Þ j 1 ;...;j K ðp; q; tÞ þ X K k¼1 W ðnÞ j 1 ;...;j k þ1;...;j K ðp; q; tÞ " # þ n " Â 0 W ðnÀ1Þ j 1 ;...;j K ðp; q; tÞ þ X K k¼1 j k k Âk W ðnÞ j 1 ;...;j k À1;...;j K ðp; q; tÞ; ð28Þ where È ¼ @=@p and These equations have been derived from factorized initial conditions 25,28,29,51) but are equally valid in the correlated case. 24)In the last reference a method to reduce the number of Matsubara frequencies needed in the numerics has been detailed, which is not employed here, due to the relatively high temperature (300 K) that will be considered.The above equations are then truncated by using the "terminators" expressed in the Wigner representation.Note that a discussion of the terminator in the density matrix case, together with a graphical representation in terms of K-faces of K þ 1 simplexes is given in Ref. 52.The number of Matsubara frequencies to be included in the calculation, K, is chosen to satisfy K ) ! c = 1 , with !c being a characteristic frequency of the system.In case that the quantity we truncate the hierarchy equations by replacing Eq. ( 28) with @ @t W ðnÞ j 1 ;...;j K ðp; q; tÞ ¼ Àð LQM þ Ä0 ÞW ðnÞ j 1 ;...;j K ðp; q; tÞ: We can evaluate W ðnÞ j 1 ;...;j K ðp; q; tÞ through numerical integration of the above equations.only the first element Wðp; q; tÞ W ð0Þ 0;0;...;0 ðp; q; tÞ has a physical meaning and the other elements W ðnÞ j 1 ;...;j K ðp; q; tÞ are initially introduced to avoid the explicit treatment of the inherent memory effects, it has been shown, however, that these elements allow us to take into account the system-bath coherence, 22) entanglement [53][54][55] and expectation values that include the bath operators as h ĤI i Àh q P j xj i. 23) The HEOM consist of an infinite number of equations, but they can be evaluated with the desired accuracy by depicting the asymptotic behavior of the hierarchal elements for different N and using this to determine whether or not there are sufficiently many members in the hierarchy.Essentially, the error introduced by the truncation turns out to be negligibly small if N is sufficiently large, which may be the case even for values lower than indicated in the inequality (32).

Numerical details
The HEOM have been studied numerically using a finite mesh representation of the Wigner function.The number of grid points in the q and p direction for the double barrier system described above are 176 and 200, respectively.For the spatial derivative in the kinetic term of the Liouvillian, a third-order left-or right-handed (depending on the sign of the momentum) upwind differencing scheme and for the second derivatives with respect to p a fourth order centered difference scheme is appropriate. 29)Simultaneously with the Wigner function, we determine the self-consistent part of the potential U self ðq; tÞ ¼ ÀeÈðq; tÞ by solving the Poisson equation, Eq. ( 2), with n À ðq; tÞ ¼ R dp Wðp; q; tÞ and the given doping density.Furthermore, the inflow boundary conditions are determined from a HEOM propagation of a free particle with periodic boundary conditions using a canonical distribution as a temporal initial state for obtaining the equilibrium state at the boundaries (see Fig. 4) of the free particle under coupling to the heat bath.
The time-step for the integration of the differential equations is chosen as 8:27 Â 10 À2 fs and we used an explicit fourth order Runge-Kutta method for time integration.To obtain the asymptotic averaged current (see below), propagation was done for a total time span of around 2000 fs.Furthermore, we fixed the frequency of the periodic gate fields to be !¼ 0:01 eV=ħ (% 1:5 Â 10 13 rad=s).
Finally, we verified our numerical results to be presented below by running calculations with different combinations of the hierarchy numbers N ¼ ð2; 3; 4Þ and Matsubara frequencies K ¼ ð2; 3Þ for the hierarchy termination and found that the (converged) results for the combination (N ¼ 4, K ¼ 3), leading to 69 additional hierarchy equations (in addition to the equation for W ð0Þ 0;0;...;0 ), did not deviate by more than approximately 1 percent from that with 9 hierarchy equations (N ¼ 2, K ¼ 2).

Electron Pumping in the Time-Domain in Presence of
Coupling to a Heat Bath We now present numerical results for the current as well as the average current through the electron pumping device in presence of a dissipative environment.We carried out the calculations using the self-consistently determined potential and, for reasons of comparison, also without self-consistent potential.

Current as a function of time with and without U self
As known from the double delta-barrier model case, the electron pumping effect depends on the phase shift φ between the two gate fields.Therefore, we first present the timedependent current divided by unit area at a temperature of 300 K, for different values of the phase difference in Fig. 5 with a gate-field amplitude of U g ¼ 0:1 eV and dissipation parameters =2 ¼ 72:5 GHz, =2 ¼ 24:2 THz.There is a vanishingly small current for ' ¼ 0 and a maximum amplitude current is observed for ' ¼ =2.
To highlight the importance of the self-consistent treatment, we also plot the current as a function of time for a gatefield amplitude of U g ¼ 0:1 eV and dissipation parameters =2 ¼ 72:5 GHz, =2 ¼ 24:2 THz with and without the self-consistent field in Fig. 6 in the case of phase difference =2.The effect of self-consistency is reducing the maximum value (and also the average value, see below) by a factor of 2!
To shed some more light on this fact we take the temporal average of the total potential defined in Eq. (1) over a (large) integer number of periods 2=! of the external forcing (see Fig. 7).The result is compared with the temporal average of  the potential without the self-consistent term, which is merely the static potential displayed in panel (b) of Fig. 1.The tunneling current is suppressed in the self-consistent case, which adds additional height to the total average barrier, making tunneling less probable than in the case without the self-consistent field.

Asymptotic average current
From the numerical calculations of the time-dependent current, e.g., displayed in Fig. 5, the asymptotic average current is obtained by taking the average over an integer number of periods of the observed current oscillations in the asymptotic regime after a time about 200 fs.To study the influence of a dissipative environment on the pumping, we have performed the average current calculations for the same gate-field parameters and the same temperature as in Fig. 5 as a function of the phase difference in the interval ' 2 ½0; and for several different values of coupling strength ζ for a fixed value of the Drude cutoff parameter =2 ¼ 24:2 THz.The corresponding numerical results are plotted in Fig. 8.There we only plot the φ-range form zero to π, because the corresponding extended curve from 0 to 2 is inversion symmetric around the point at π (see Fig. 3).It is interesting to note that in the time-domain, we observe (not shown) that the current at ' ¼ 0 is exactly zero (apart from numerical noise), while the current at ' ¼ is nonzero but averages to zero.
By increasing the dissipation strength, we observe a decrease of the average pumped current, leading to a dissipation-induced decrease in the maximum pumped current through the resonant tunneling structure.We note that there are regions in the parameter space where an increase of dissipation can also lead to an increase in the maximum tunneling current (not shown).The behavior of the maximum as a function of the bath coupling (and Drude parameter) is non-monotonic.In contrast to the dissipationinduced enhancement at ' ¼ 0, reported in Fig. 3(b) of Ref. 20, here, the maximum value of the current near ' ¼ =2 is influenced by the dissipation.In Ref. 20 dissipation has led to the shift from a sine-like curve to a cosine-like curve as a function of φ without(!) a change of the maximum current.Furthermore, it was shown for the case of a periodic potential without self-consistent force that the ratchet current is a decreasing function of the dissipation, while the value of the phase difference for the maximum current approaches towards =2. 25)he calculated average current is also greatly modified with or without self-consistent field, as displayed in Fig. 8(b).Firstly, it is larger by a factor of 2 compared to the selfconsistent case and secondly, the position of the maximum value (as a function of φ) as well as the maximum value itself does not change appreciably as a function of the damping strength in the case without self-consistent potential.

Conclusions and Outlook
We have studied numerically in the time-domain the buildup of an average electronic dc current in a double barrier quantum well structure through the application of symmetry breaking external gate fields.To this end, we have used an approximation-free non-Markovian and non-perturbative propagation technique of the Wigner function in phase space.
Taking into account the coupling of the electron dynamics to a reservoir of phonon modes with a Drude-Lorentz  spectral density in the framework of reduced we were able to show that the maximum current that is generated by the gate fields (near a phase difference of =2) is changing (non-monotonically) as the coupling strength to the phonon modes increases.The overall effect of dissipation is stronger in the case of self-consistent determination of the potential.Furthermore, also qualitatively, self-consistency leads to more asymmetric behavior of the current near ' ¼ =2 in comparison to the non-self-consistent calculation.This finding extends recent results for different physical systems (2 site molecular wires under bi-harmonic driving 20) and Cooper pair sluice 21) ).In contrast to the results found in Ref. 20, in which a mere shift of the current versus phase difference was observed, here we see a change of the maximal pumped current with increasing dissipation.
Using the HEOM methodology, it would be worthwhile to investigate also the case of more than two barriers, leading to several quantum wells and several phase differences that can be varied to generate a net dc current.A recent study in the dissipation-less case 38) illustrates interesting reversals of the pumped current's direction, and it would be interesting to elucidate how this behavior is influenced by phonon coupling.Furthermore, a much more demanding but nevertheless possibly very fruitful direction of future research would be the calculation of correlated transport of electrons through quantum dots.Finally, the calculation of heat currents is closely related to the electron current in the adiabatic case. 10)It may also be calculated in the timedependent fashion that we have used here.

Fig. 1 .
Fig. 1. (Color online) (a) Resonant tunneling structure with the quantum well of undoped GaAs and two barriers of undoped AlGaAs; spacer layers consist of undoped GaAs and contact regions of doped GaAs (doping concentration of 2 Â 10 18 cm −3 ); (b) static double barrier potential.

Fig. 4 .
Fig. 4. (Color online) Representation of the inflow boundary conditions for electron transport through a device.