A nonequilibrium theory for transient transport dynamics in nanostructures via the Feynman-Vernon influence functional approach

In this paper, we develop a nonequilibrium theory for transient electron transport dynamics in nanostructures based on the Feynman-Vernon influence functional approach. We extend our previous work on the exact master equation describing the non-Markovian electron dynamics in the double dot [Phys. Rev. B78, 235311 (2008)] to the nanostructures in which the energy levels of the central region, the couplings to the leads and the external biases applied to leads are all time-dependent. We then derive nonperturbatively the exact transient current in terms of the reduced density matrix within the same framework. This provides an exact non-linear response theory for quantum transport processes with back-reaction effect from the contacts, including the non-Markovian quantum relaxation and dephasing, being fully taken into account. The nonequilibrium steady-state transport theory based on the Schwinger-Keldysh nonequilibrium Green function technique can be recovered as a long time limit. For a simple application, we present the analytical and numerical results of transient dynamics for the resonance tunneling nanoscale device with a Lorentzian-type spectral density and ac bias voltages, where the non-Markovian memory structure and non-linear response to the bias voltages in transport processes are demonstrated.


I. INTRODUCTION
The investigation of quantum coherence dynamics far away from equilibrium in nanoscale electronic devices has been attracting much attentions in the past decade due to various applications in nanotechnology and quantum information processing. Experimentally it became possible not only to directly manufacture structures but also to investigate their nonequilibrium quantum coherence properties under well controlled parameters. [1][2][3][4][5][6][7] Theoretically, tremendous studies have been done focusing on the time-independent steady-state to analyze the nonequilibrium transport properties. However, an increasing number of work employing various theoretical techniques (see Refs. [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24]) has recently occurred addressing the timedependent electron dynamics far away from the equilibrium. In fact, for practical applications, a full understanding of nonequilibrium dynamics under external field controls, i.e., the time evolution of electrons in nanoelectronics devices from some initial preparation toward any specifically designed state within an extremely short time, has been receiving more and more attentions.
The prototypical nanostructure concerned in this paper is made by a gate-defined region on a semiconductor containing a few of discrete electronic states. Electrodes are implanted around this central region to control the electrons across it. Also gates are deposited to adjust the electronic states within the central area as well as their couplings to the surrounding electrodes. Most of the theoretical developments to nanoelectronics devices have been focusing on the understanding and prediction of transport dynamics in these devices. However, as a quantum device, in particular for quantum information processing, 4-7 the big challenge is to understand and predict not only how fast or slow a nanoelectronic device can turn on or off a current, but also how reliable and efficient the device can manipulate quantum coherence of the electron states through external bias and gate voltage controls. This requires a nonequilibrium quantum device theory to analyze the time-dependent quantum transport intimately entangling with quantum coherence of electrons inside the device. It is the purpose of this paper to attempt to establish a nonequilibrium theory of transient dynamics for electron quantum states accompanied with electron quantum transport in a nanoelectronic device that can be feedbackly controlled through the nonlinear response to the external time-dependent bias and gate voltage pulses.
Physically the nanostructure we are going to study is typically an open quantum system in the sense that its central region has exchanges of particles, energy and information with its surroundings. From an open quantum system point of view, the nonequilibrium electron dynamics is completely described by the master equation through the physical observable O(t) = tr[Oρ(t)], where ρ(t) is the reduced density matrix of the central system that can fully depict the dynamics of electron quantum states and is determined by the master equation. The transient electron transport should be in principle solved from the reduced density matrix as well. However, it has been struggled for many years without a very satisfactory answer in finding the exact master equation for an arbitrary open quantum system. 25 Most of the master equations used in the literature are obtained using semiclassical approximation or perturbation truncation, such as the semiclassical Boltzmann equation 26,27 or master equations under the Born-Markov approximation. 28,29 However, for the nanostrucutres with the extremely short length scales (∼1 nm) and extremely fast time scales (∼1 fs), the semiclassical Boltzmann equation or the master equation under the Born-Markov approximation are most unlikely applicable. 24 Literarily, there are two fundamental but equivalent methods in dealing with modern nonequilibrium physics. They are the Schwinger-Keldysh nonequilibrium Green function technique [30][31][32] and the Feynman-Vernon influence functional approach. 33 In the practical applications, both methods have their own advantages. The Schwinger-Keldysh nonequilibrium Green function technique is a standard method for systematic perturbation calculations to various nonequilibrium physical processes in many-body systems. 34,35 As it is evidenced the nonequilibrium Green function technique provides an extremely powerful tool in the study of the steadystate transport properties in nanostructures. 24,[36][37][38] In the steady-state limit where the initial state of the nanostructure is irrelevant, the transport current can be expressed in terms of the nonequilibrium Green functions in the frequency domain which largely simplifies the transport problem. Without the time-dependence, the physics of the nonequilibrium transport is essentially determined by the density of states (the difference between the retarded and the advanced Green functions). The nonequilibrium Green function technique has been used to successfully describe a variety of new phenomena, such as Kondo effect, Fano resonance and Coulomb blockade effects in quantum dots. The extension of the nonequilibrium Green function techniques to the time-dependent transport phenomena in nanostructures has also been developed [10][11][12][13]24 where the time-dependent nonequilibrium Green functions must be maintained. Except for the wide band limit (WBL), 11,24 the relevant Green function calculations become rather complicated in order to taking into account the different initial quantum state effect and the non-Markovian memory structure. 17,20 On the other hand, the Feynman-Vernon influence functional approach 33 are mainly used to study dissipation dynamics in quantum tunneling problems 39 and decoherence problems in quantum measurement theory. 40 It is very powerful to derive nonperturbatively the master equation for the reduced density matrix of an open system by integrating out completely the environmental degrees of freedom through the path integral approach, where the non-Markovian memory structure is manifested explicitly. In the early applications of the influence functional approach, the master equation was derived for some particular class of Ohmic (white-noise) environment for the quantum Brownian motion (QBM) modeled as a central harmonic oscillator linearly coupled to a set of harmonic oscillators simulating the thermal bath. 41,42 The exact master equation for the QBM with a general spectral density (color-noise environments) at arbitrary temperature that can fully address the non-Markovian dynamics was the Hu-Paz-Zhang master equation. 43 Applications of the QBM exact master equation cover various topics, such as quantum decoherence, quantum-toclassical transition, quantum measurement theory, and quantum gravity and quantum cosmology, etc. 29,44,45 Very recently, such an exact master equation is further extended to the systems of two entangled bosonic fields [46][47][48] for the study of non-Markovian entanglement dynamics. 49 Nevertheless, using the influence functional approach to obtain the exact master equation has been largely focused on the bosonic (thermal) environments in the past half century.
Meanwhile, using master equation to study quantum electron transport has also been attracting much attentions recently. With the nonequilibrium Green function technique, the master equation can be formally obtained in terms of real-time diagrammatic expansion up to all the orders, 9,[50][51][52] although practically most of the master equation approach used in quantum transport by far only take the perturbative theory up to the second order. 15,[53][54][55][56][57] On the other hand, with the help of the influence functional approach, an alternative formalism for the master equation for studying quantum transport and electron dissipative dynamics has recently been developed in term of the hierarchical expansion. 18 This new master equation approach provides a very efficient tool for the numerical study of the quantum transport properties, including the accurate evaluations of the Coulomb blockage and Kondo transition dynamics. 23,58 In this paper, we shall develop a nonequilibrium theory to describe nonperturbatively the transient electron transport dynamics intimately entangling with the electron quantum decoherence in nanostructures. By extending the influence functional approach to the fermionic reservoirs, we have recently derived an exact master equation for a large class of nanostructures. 21 This exact master equation is capable for the study of the full non-Markovian decoherence dynamics of electrons in nanoelectronics devices. In the present work, we shall extend the previous work to the nanostructures in which the energy levels of the central region, the couplings to the leads and the external biases applied to leads are all timedependent. We then derive the exact transient current in terms of the reduced density matrix within the same framework. The resulting transient current provides an exact non-linear response theory that can be used to investigate the transient quantum transport processes accompanied explicitly with the non-Markovian quantum relaxation and dephasing dynamics. The nonequilibrium steady-state transport theory based on the nonequilibrium Green function technique that is widely used in the literature can be recovered as a long time limit from the present theory.
The remainder of this paper is organized as follows.
In Sec. II, we outline the derivation of the exact master equation we obtained recently, 21 with the extension to the time-dependent Hamiltonian for a general nanostructure. In Sec. III, we derive in detail the exact transient current in terms of the reduced density operator based on the Feynman-Vernon influence functional approach. The resulting transient current I α (t) through the lead α is given by the general formula, Eq. (19), or equivalently Eq. (30), which is valid for arbitrary bias and gate voltage pulses with arbitrary energy-dependent couplings between the central system and the leads. The relation between the reduced density matrix and the transient current is established through the master equation, see Eq. (23) where the superoperators L + α (t) and L − α (t) given explicitly by Eq. (20) and (22) encompass all the backreaction effects associated with the non-Markovian dynamics of the central system interacting with the lead α.
The steady-state current formula based on the nonequilibrium Green function technique can be recovered as a long time limit, see Sec. III.E. Applications to transient dynamics with Lorentzian-type spectral density and various time-dependent ac bias voltages are demonstrated in Sec. IV. The summary and prospective are given in Sec. V. In Appendix, a close connection between the Feynman-Vernon influence functional approach and the Schwinger-Keldysh nonequilibrium Green function technique for quantum transport theory is established.

II. EXACT MASTER EQUATION FOR NANOSCALE ELECTRONIC DEVICES
In this section, we shall outline the exact master equation derived recently by two of us for nanoscale electronic devices 21 with the extension to the time-dependent energy levels and couplings. We shall only list the necessary formula that we will use later for the derivation of the transient current. For more detailed derivation, please refer to Ref. [21]. We begin with the Hamiltonian of a prototypical mesoscopic system for electron transport, Here the first summation is the electron Hamiltonian H S for the central system where the electron-electron interaction is not considered. But as we have shown we can properly and explicitly exclude the doubly occupied states and the resulting master equation can be applied to the strong Coulomb repulsion (i.e. Coulomb blockage) regime. 21 The second summation is the Hamiltonian α H α describing the noninteracting electron leads (the source and drain electrodes, · · · ) labeled by the index α (= L, R, · · · ). The last term is the electron tunnel-ing Hamiltonian H T between the leads and the central system. Since we will focus on the transient dynamics of quantum transport in this paper, the bias voltages V α applied to leads are considered to be time-dependent. Thus the single electron energy levels in the leads is changed as ǫ αk → ǫ αk (t) = ǫ αk + eV α (t). Meanwhile, the energy levels ǫ ij of the central system and the couplings V iαk between the central system and the leads are controllable through the gate voltages and external fields so that they can be, in general, also time-dependent: Throughout this work, we set = 1, except for the transient current where we will put back from its definition.
The master equation for the nonequilibrium electron dynamics of the central system is derived based on the Feynman-Vernon influence functional approach. 33 As it is well known, 25 the nonequilibrium electron dynamics of an open system is completely determined by the socalled reduced density matrix. The reduced density matrix is defined from the density matrix of the total system (the central system plus the leads, and the leads are treated as the reservoirs to the central system) by tracing over entirely the environmental degrees of freedom: ρ(t) ≡ tr R ρ tot (t), where the total density matrix is formally given by ρ and T is the time-ordering operator. As usual, we assume that the central system is uncorrelated with the reservoirs before the tunneling couplings are turned on: 39 ρ tot (t 0 ) = ρ(t 0 ) ⊗ ρ R (t 0 ), and the reservoirs are initially in the equilibrium state: ρ R (t 0 ) = 1 Z e − P α βα(Hα−µαNα) where β α = 1/(k B T α ) is the initial inverse temperature and N α = k c † αk c αk the particle number operator for the lead α. Then the reduced density matrix at an arbitrary later time t can be expressed in terms of the coherent state representation 59,60 as with ξ = (ξ 1 , ξ 2 , · · · ) andξ = (ξ * 1 , ξ * 2 , · · · ) being the Grassmannian numbers and their complex conjugate defined through the fermion coherent states: a i |ξ = ξ i |ξ and ξ|a † i = ξ|ξ * i . The propagating function in Eq. (2) is given in terms of Grassmannian number path integrals, where S c [ξ, ξ] is the action of the central system in the fermion coherent state representation, and F [ξξ;ξ ′ ξ ′ ] is the influence functional obtained after integrating out all the environmental (reservoirs) degrees of freedom: The two-time correlations in Eq. (4) are defined by which depict all the time-correlations of electrons in the leads through the interactions with the central system, and f α (ǫ αk ) = 1/(e βα(ǫ αk −µα) − 1) is the Fermi distribution function of the lead α at the initial time t 0 . This influence functional takes fully into account the back-reaction effects of the reservoirs to the central system. It modifies the original action of the system into an effective one, e which dramatically changes the dynamics of the central system. The detailed change is manifested through the generating function of Eq. (3) by carrying out the path integral with respect to the effective action S eff [ξξ;ξ ′ ξ ′ ]. While the path integral D[ξξ;ξ ′ ξ ′ ] integrates over all the forward pathsξ(τ ), ξ(τ ) and the backward pathsξ ′ (τ ), ξ ′ (τ ) in the Grassmannian space bounded byξ(t) =ξ f , ξ(t 0 ) = ξ 0 andξ ′ (t 0 ) =ξ ′ 0 , ξ ′ (t) = ξ ′ f , respectively. Since S eff [ξξ;ξ ′ ξ ′ ] is only a quadratic function in terms of the integral variables, the path integrals in Eq. (3) can be reduced to Gaussian integrals so that we can use the stationary path method to exactly carry out these path integrals. 61 The resulting generating function is: 21 in which the time-dependent coefficients are given explicitly as: Here, we have also expressed the stationary paths in terms of N × N matrix variables u(τ ), v β (τ ) andū(τ ) for t 0 ≤ τ ≤ t, where N is the total number of single particle energy levels in the central region, and the index β implies that the corresponding function is temperature dependent. They satisfy the following dissipationfluctuation integrodifferential equations (i.e. the station-ary path equations of motion): with the boundary conditions u(t 0 ) = 1,ū(t) = 1 and v β (t 0 ) = 0, where g α (τ, τ ′ ) and g β α (τ, τ ′ ) are the nonlocal two-time correlation matrix functions, whose matrix elements are given by Eq. (5). As we will show in Appendix the matrix variables u(τ ),ū(τ ) and v β (τ ) correspond respectively to the retarded, advanced and lesser Green functions, and g α (τ, τ ′ ) and g β α (τ, τ ′ ) are the retarded and lesser self-energies in the nonequilibrium Green function technique. For the most general cases where all the parameters in the Hamiltonian (1) are time-dependent, the time translational invariance is broken down. Then u(τ ) andū(τ ) are usually independent except for the end-points where we haveū(t 0 ) = u † (t). If all the parameters in the Hamiltonian of Eq. (1) are the time-independent, u(τ ) will become only a function of τ − t 0 andū(τ ) = u † (t − τ + t 0 ), as we have shown in Ref. [21].
Taking the time derivative of the reduced density matrix with the solution of the propagating function (6), together with the D-algebra of fermion creation and annihilation operators in the fermion coherent state representation, we arrive at the final form of the exact master equation we obtained previously 21 The first term (the commutator) in the master equation accounts for the renormalized effect (including the timedependent shifts of the energy levels and the changes of the transition amplitudes between them) of the central system due to the interaction with the leads. The resulting renormalized Hamiltonian is H ′ While the rest terms in the master equation describe the dissipation and noise effects (which results in a non-unitary evolution of the central system) induced also by the interaction with the leads. All the timedependent coefficients in Eq. (8) are determined explicitly by u(t), u † (t) and v β (t) as follows: Here we have used the relations obtained from Eq. (7), As it is well known, the time-dependent coefficients in the master equation describe the non-Markovian dynamics due to the interaction between the central system and the leads. The time-dependence of these coefficients is fully determined by the functions u(τ ),ū(τ ) and v β (τ ) which are governed by the integrodifferential equations (7) where the time integrals involve the nonlocal time correlation functions of the reservoirs, g α (τ, τ ′ ) and g β α (τ, τ ′ ). These non-local time correlation functions characterize all the non-Markovian memory structures of the central system interacting with its environment through the coupling Hamiltonian H T . By solving Eq. (7), one can completely describe the quantum decoherence dynamics of electrons in the central region due to the entanglement between the central system and the leads. Thus, the master equation determines indeed the exact nonequilibrium dynamics of the system. It may be also worth mentioning that our master equation does not need to be in Lindblad form to preserve positivity since it is exact so that the dampening coefficients are time-dependent and the positivity is guaranteed. In the next section, we will derive the transient current passing through the central system within the same framework. The result will allow us to explore the transient dynamics of the nonequilibrium quantum transport accompanied with the non-Markovian decoherence phenomena, which receives more and more attentions recently due to rapid developments in nano-technology, spintronics and quantum information processing, etc. 6,7

III. EXACT TRANSIENT CURRENT
A. Influence functional approach The transient current from the α-lead tunneling through α-junction into the central region is defined in the Heisenberg picture as Here, e is the electron charge, , and ρ H tot is the total density matrix in the Heisenberg picture. By explicitly calculating the above commutation relation with the Hamiltonian of Eq. (1) and then transforming it into the Schördinger picture, we have where the operators The time-dependence of these two operators comes from the non-Markovian memory dynamics by tracing over the environmental degrees of freedom of the total density matrix ρ tot (t) in the Schrödinger picture. These two operators are indeed the effective (dressed) electronic creation and annihilation operators acting on the reduced density matrix of the central system, as a result of tracing over the reservoir degrees of freedom for the corresponding operators acting on the lead α.
Following the procedure of obtaining the reduced density matrix through the influence functional approach, 21 we can write where the operator-associated propagating function is defined as Similar to the calculation of the influence functional in Eq. (14), the operator-associated influence functional F A αi [ξξ;ξ ′ ξ ′ ] can be calculated in the same way with the result, where the functional expression of the effective electron annihilation operator is obtained as and the same influence functional F [ξξ;ξ ′ ξ ′ ] is given by Eq. (4). It should be pointed out that the factorability of the operator-associated influence functional is due to the fact that the environmental path integral is exactly computable.
Since the path integrals in Eq. (14) can also be reduced to Gaussian integrals, similar to the derivation of the master equation for the reduced density matrix (8), we use again the stationary phase method to exactly carry out the path integrals of Eq. (14). The result is: where is given by Eq. (6) in the last section. Substituting Eq. (17) into Eq. (13), and using the identity, together with the D-algebra for the fermion creation and annihilation operators in the fermion coherent state representation, we obtain the effective electronic annihilation operator after tracing over completely the environmental degrees of freedom: where the coefficient matrices λ β α (t) and κ α (t) are the same as these appeared in the master equation (8) and are explicitly given by Eq. (10), ρ(t) is just the reduced density matrix determined by the master equation.
Accordingly, substituting the above result into Eq. (12), the transient current can be directly calculated. The resulting transient current is rather simple: where the notation Tr is the trace over the energy level basis of the central system. ρ )] is the single particle reduced density matrix. While the matrix elements of g α (τ, τ ′ ) and g β α (τ, τ ′ ) are defined by Eq. (5), and u(τ ), v β (τ ) andū(τ ) are determined nonperturbatively by Eq. (7), as shown in the last section. Thus the transient current I α (t) that flows from the lead α is completely determined within the framework of the master equation for the reduced density matrix. As we shall show later, from the above result we can easily reproduce the steady-state current at t → ∞ in terms of the nonequilibrium Green functions and the Landauer-Büttiker formula.

B. Relation between the transient current and the reduced density matrix
In fact, the transient current defined by Eq. (12) can be written as a trace over the transient current operator: , where the current operator is obtained directly from Eq. (18): (20) Here we have introduced a current superoperator L + α (t) such that the current operator can be simply expressed as the current superoperator acting on the reduced density matrix. Note from Eq. (10) that α (λ β α +λ β † α ) = γ β and α (κ α + κ † α ) = 2γ, the transient current operator given by Eq. (20) is closely connected to the master equation (8) for the reduced density matrix ρ(t).
Because of the trace over the states of the central system, the transient current of Eq. (12) can be alternatively expressed as Then using Eq. (18) again, we havễ where L − α (t) is defined as the superoperator forÎ α (t). Now it is easy to check that the master equation (8) for the reduced density matrix and the transient current of (19) can be simply expressed as where H S is the original Hamiltonian of the central system. This analytical operator relation between the reduced density matrix and the transient current shows explicitly the intimate connection between quantum decoherence and quantum transport in nonequilibrium dynamics. The reservoir-induced non-Markovian relaxation and dephasing in the transport processes are manifested through the superoperators of (20) and (22) acting on the reduced density matrix ρ(t). The time-dependent parameters λ β α (t) and κ α (t) appeared in the superoperators are given by Eq. (10) which are determined by the dissipation-fluctuation integrodifferential equation (7). This completes the nonequilibrium theory for transient electron dynamics in nanostructures we concerned. It should be pointed out that Eq. (23) was formally expressed in terms of the real-time diagrammatic expansion and the hierarchical expansion. 9,18 Here we have the analytical solution.

C. Current conservation law
Now we should derive the current conservation law for a self-consistent check. Consider the equation of motion for the operator a † i a j in the Heisenberg picture, Taking the expectation value of the above equation with respect to the state ρ H tot (the total density matrix in the Heisenberg picture) and then transforming it into the Schördinger picture, we obtain Here we have used the definition of the single-particle reduced density matrix again, ρ , and also introduced a current matrix Equivalently, the transient current of Eq. (12) is simply given by I α (t) = e TrI α (t). In other words, the equation of motion for On the other hand, the equation of motion for the single particle reduced density matrix can also be obtained easily from the exact master equation (8). The result is Using the expression for γ β (t) given by (9c), and comparing Eqs. (25) with (26), we have The above equation reproduces exactly the transient current of Eq. (19). This also provides a self-consistent check for the expression of transient current derived directly from the influence functional approach. Taking the trace over the both sides of Eq. (25) and also noting the fact that N (t) = Trρ (1) (t), the total electron occupation in the central system, we have where I dis (t) is defined as the transient displacement current. 62 It tells that sum over the currents flowing from all the leads into the central region equals to the change of the electron occupation in the central region, as we expected. In the steady-state limit t → ∞, N (t) = Trρ (1) (t) = 0 so that I dis (t) = 0, as a consequence of current conservation.
D. Solution to single particle reduced density matrix and transient current Furthermore, we can rewrite Eq. (9c) as Comparing between Eq. (25) for the single particle reduced density matrix ρ (1) (t) and the above equation for v β (t), it indicates that the solution of ρ (1) (t) is just v β (t), apart from an initial function. Explicitly, Eq. (25) and the above equation lead to It is easy to find from the above equation the following relationship between ρ (1) (t) and v β (t): where ρ (1) (t 0 ) is the initial single particle reduced density matrix. Accordingly, the transient current (19) is reduced to The last term shows the explicit dependence on the initial single particle reduced density matrix (including the initial electron occupation in each level and the electron quantum coherence between different levels in the central region). This term is an important ingredient in the study of transient dynamics for practical manipulation of a quantum device in real time and it will usually vanish in the steady-state limit t → ∞.

E. Steady-state limit
To make an explicit comparison with the steady-state current in terms of the non-equilibrium Green functions and the Landauer-Büttiker formula used in the literature, we introduce the spectral density of the lead α: Γ αij (ω) = 2π k V αik V * αjk δ(ω − ǫ αk ) and take a timeindependent bias voltage explicitly. Then the two-time correlation functions of the α-lead can be written as where f α (ω) = 1/(e βα(ω−µα) +1) is the Fermi distribution function of the α-lead at the initial time t 0 . Using the Laplace transformation, i.e., i.e., g α (ω) is the retarded self-energy induced from the lead α. The Laplace transformation of Eq. (7a) for u(t) gives where G r (ω) is just the retarded Green function, and Σ r (ω) sums over Σ r α (ω) for all α. The advanced Green function is simply given byū(ω) = −iG a (ω) = u † (ω). Furthermore, for the time-independent Hamiltonian, the explicit solution of Eq. (7c) is Taking the long time limit: t → ∞, its Laplace transformation gives Here we have also used the relation g β (ω) = −iΣ < (ω). More explicit relations between u(t),ū(t) and v β (t) with the retarded, advanced and lesser Green functions in the real-time domain are presented in Appendix.
Substituting above results into Eq. (30) in the steadystate limit t → ∞, we obtain the steady-state single particle reduced density matrix and current: This reproduces the steady-state current in terms of the nonequilibrium Green functions in the frequency domain that has been widely used. If we consider specifically a system coupled with left (source) and right (drain) electrodes, i.e. α = L and R, respectively, and also assume that the spectral densities for the left and right leads have the same energy dependence: Γ L (ω) = λΓ R (ω), where λ is a constant. Then the net steady-state current flowing from the left to the right lead is given by This is the generalized Landauer-Büttiker formula. 24 Thus, we have simply recovered the nonequilibrium transport theory at the steady-state limit. In fact, we can also reproduce the transient transport theory in terms of the nonequilibrium Green function technique, as shown in Appendix. We now summarize the main results derived in this Section. We obtain the general formula of the transient current I α (t) through the lead α which is given by Eq. (19) or equivalently Eq. (30), accompanied with the exact master equation for the reduced density matrix ρ(t). This general transient current is valid for arbitrary bias and gate voltage pulses with arbitrary couplings between the central system and the leads. We also establish explicitly the connection of the transient current with the reduced density matrix, i.e. Eq. (23), through the superoperators L + α (t) and L − α (t) determined by Eq. (20) and (22), which encompass all the backreaction effects associated with the non-Markovian dynamics of the central system interacting with the lead α. The current conservation law is also explicitly derived. The steady-state current formula based on the nonequilibrium Green function technique is reproduced as a long time limit. These results allow us to explicitly analyze the time-dependent quantum transport phenomena intimately entangling with the electron quantum coherence and non-Markovian dynamics through the reduced density matrix. The latter describes completely the evolution of electron quantum coherence inside the nanoelectronics device. Therefore, once one solves the equation of motion (7) for u(τ ),ū(τ ) and v β (τ ), the full nonequilibrium dynamics of the nanostructures can be analyzed explicitly.

IV. ANALYTICAL AND NUMERICAL ILLUSTRATIONS
A. Time-independent bias voltage in the WBL: an analytic solution for transient dynamics For practical applications, we take a simple system as an illustration: a single dot containing only one spinless level that couples to the left and right leads. The Hamiltonian of the dot is simply written as H = ǫa † a. This system only contains two states, the empty state |0 and the occupied state |1 . All the corresponding matrixes (denoted by the bold symbols) in the above formula are then reduced to a single function, such as u(t) = u(t), v β (t) = v β (t), and ρ (1) (t) = 1|ρ|1 = ρ 11 (t) = N (t). We will first consider a time-independent bias voltage V . For simplicity, we also assume that the tunneling couplings between the leads and the dot as well as the densities of states for the leads are energy-independent. In other words, the spectral density becomes a constant Γ α (ω) = Γ α . The non-local time correlation functions is reduced to This corresponds to the wide band limit (WBL) in the literature.
To be explicit, we take the initial time t 0 = 0 and let e = = 1 in the following calculations. In the WBL, the solution of Eq. (7) is where Γ = Γ L + Γ R and v β st is the solution of v β (t) at the steady-state limit: st .
The electron occupation of the dot is calculated by Eq. (29) as where N (0) = ρ (1) (0) is the initial electron occupation of in the dot. Obviously, in the steady limit, N st = v β st . This is not surprised since ρ (1) (t) and v β (t) obey the same equation of motion that must lead to the same result in the steady-state limit.
The transient current can then be analytically obtained: where the steady-state current is Due to charge conservation, it is necessary to check the displacement current which obeys the relation (28). In the WBL, At the steady-state limit, the displacement current I dis (t → ∞) = 0, as a consequence of current conservation. It is also straightforward to calculate the net current: where the stationary net current is As we can see, once we solve the dissipation-fluctuation integrodifferential equations, Eq. (7), a complete timedependence of all the physical quantities, such as the electron occupations in the central system and the currents flowing from each lead to the central region, can be obtained with the explicit dependence of the time and the initial electron occupation without ambiguity. From Eq. (42), we see that the initial current I α (0) = −Γ α N (0) which depends on the initial occupation of the dot. This result is consistent with the electron occupation in the dot, Eq. (41). For zero initial occupation, the initial current is zero. Note that some of the above results are also obtained recently using the nonequilibrium Green function technique 20 .

B. Non-Markovian memory structure
Realistically, the spectral density of the leads must depend on the energy. Here we take the energy dependence as a Lorentzian-type form: 17,18,21 where Γ α describes the coupling strength and W α is the line width of the source (drain) reservoir with α = L(R).
Obviously the WBL, Γ α (ω) = Γ α , is achieved by simply letting W α → ∞. The lead correlation functions with time-independent voltage can be parameterized as 18 The first term in Eq. (48b) with m = 0 arises from the pole of the spectral density function, with The other terms with m > 0 (M → ∞ in principle) arise from the Matsubara poles, where the relevant parameters are explicitly given as There are four typical timescales in an open system to characterize the non-Markovian memory structure: the timescale of the central system (∼ 1/ǫ), the timescale of the reservoirs (∼ 1/W α ), the mutual timescale due to the coupling between the system and the reservoir (∼ 1/Γ α ), and the thermal timescale (∼ β α or 1/µ α ). The timescale of the central system is usually fixed when the system is set up. Then comparing with the character time of the central system, a smaller finite line width W α , a relatively large coupling Γ α , a low temperature 1/β α , and a comparable bias voltage eV = µ L − µ R to the energy levels of the central system will lead to a stronger non-Markovian memory processes for a Lorentzian-type spectral density, as we have demonstrated in [21]. In fact, the line width W α in a Lorentizan-type spectral density is the main factor leading to the non-Markovian dynamics in transient transport. In the WBL, W α → ∞, the dominated memory structure are mostly washed out. This can be seen directly from the reservoir correlation functions. The correlation function g α (t − τ ) of Eq. (48a) can be simplified to Γα 2 δ(t − τ ) in the WBL. While for g β α (t − τ ) in Eq. (48b), the first term (m = 0) is also simplified to a delta function of t − τ but the other terms (m ≥ 1) are apparently changed not too much: The profile of this temperature-dependent time correlation function is plotted in Fig. 1. If we take further a high temperature limit β α → 0, the summation term in Eq. (51) will also be reduced to a delta function of t − τ , see Fig. 1. Then no any memory effect remains, and a true Markov limit is reached at high temperature limit. On the other hand, for a large bias voltage limit This reduces the electron dynamics into a Markov limit again, as it has been widely used in the literature. [63][64][65] However, for a relatively low temperature or a finite bias voltage, there appears to exist some non-Markovian effect in the WBL, coming from the summation term in Eq. (51). Fig. 2 shows the time dependence of the correlation dependence with different voltages. As one can see, the temperature-dependent time-correlation function does not approach to a delta function in time.
To examine if some non-Markovian effect can still survive in the WBL, we numerically calculate the transient current passing through the single level resonant tunneling nanostructure considered in the last subsection. In the sequential tunneling regime µ L > ǫ > µ R , the exact solutions to the occupation and the transient current are close to the Markov limit (differing by a few present except for a very short time scale from the beginning), as shown in Fig. 3(a)-(d). While, in the co-tunneling regime with µ α ≫ ǫ, the exact solution of the occupation and In other words, when the line width W → ∞, not only for the high temperature and large bias limit, but even for a finite temperature and a finite bias voltage, the non-Markovian effects becomes quite weak that are most likely negligible in experiments. Thus the WBL mainly takes into account the Markov dynamics. The manifestation of the non-Markovian memory structure then should go beyond the WBL. 21 In Fig. 4, we calculate the transient current with different line widths to demonstrate the non-Markovian effect in transport phenomena. As we can see, the transport dynamics is significantly different from the Markov limit for a small W . Increasing the value of W will decrease the memory effect accordingly. When the line width W ≥ 50Γ, the exact solution of the transient current closely approaches to the Markov result that is consistent with the WBL. A analysis of the non-Markovian dynamics in this simple resonant tunneling system has also recently been studied using Heisenberg equations of motion. 66

C. Transient transport dynamics with time-dependent bias voltage
We are now in the position to study the transient transport phenomena in response to time-dependent bias voltages. The calculation of the time-dependent electron current in response to time-dependent bias voltages has received a great attention recently, through various different theoretical approaches. 8,9,[11][12][13][14]17,19,20,23,24 Transient transport dynamics is also a central ingredient in many different experiments, such as single-electron pumps and turnstiles with time-modified gate signals moving electrons one by one through quantum dots, [67][68][69][70][71] and the study of the quantum capacitance and inductances with ac voltage response. [72][73][74][75][76] All these problems can be studied explicitly in the present theory now. The corresponding nanoscale devices can be modelled as a resonant tunneling nanostructure considered in this section. The applying time-dependent voltages are taken to be the most commonly interesting ones. For time-dependent voltages, as we mentioned in Sec. II the single particle energy levels of the leads are changed to ǫ αk (t) = ǫ αk + eV α (t). The non-local time correlation functions of the leads can be expressed as With the parametrization of Eq. (48), it is not difficult to numerically calculate the transient electron dynamics for arbitrary line width W α . In the following calculation, we set the symmetric ac voltages, i.e., µ L (t) = eV (t)/2 and µ R (t) = −eV (t)/2. For the quantum dot with a single level, the reduced density matrix can be fully characterized by the electron occupation in the dot. The exact numerical results for the transient current and occupation due to different types of applied ac voltages are presented as follows.

Exponentially time-dependent bias voltage
We shall first study the transient transport dynamics in response to an exponentially time-dependent bias voltage V (t) = V (1 − e −t/τ ), where τ > 0 is a time dominating switch-on rate of the voltage. The asymptotic limit τ → 0 + corresponds to a step function. The numerical results are plotted in Fig. 5. The first peak shown in the displacement current arises from the cotunneling process during the initially short time scale. At beginning, the Fermi surface of the both leads are nearly equivalent to zero and the dot energy level is higher than the Fermi surface. The initial currents through both leads are equal to each other due to the totally symmetric structure. This leads to a zero initial net current. This feature is common for other type of ac bias voltages discussing later. In general, the emergence of the peaks in the current corresponds to a steep transient behavior of the electron states (i.e. the occupation for the single-level dot) occurs inside the dot. The non-linear response to the time-dependent bias is clearly manifested in the change of both the electron state in the dot and the transient current through the leads (including the individual current through each lead and the displacement and net currents), as we plotted in the Fig. 5. In particular, the net current changes in time that closely follows the change of electron occupation in the dot, while the displacement current depicts the steep change rate of the occupation, as we expected from Eq. (28).  where the Bessel function satisfies J −n (z) = (−1) n J n (z) and ∆ L,R = ∓ eVc 2 .

Gaussian pulse
The last example we shall study is the transient dynamics droven by a Gaussian pulse V (t) = V exp{− (t−τ1) 2 τ 2 2 }. The width and the center of the pulse are determined by τ 2 and τ 1 , respectively. The exact numerical result plotted in Fig. 8 shows that the net current peak emerges (a slightly delay) just after the voltage pulse while the corresponding response of the electron occupation delays significantly. This behavior is easy to be understood that the external voltage pulse leads to a sharp change of the electron occupation due to the delay response. The shape change of the transient current comes from the largest change rate of the electron occupation in the dot. The delay response effect is determined by the tunneling rate Γ. Outside the voltage pulse, the tunneling current comes completely from the cotunneling effect and the occupation in the dot decays to a stationary value. It is interesting to note that the response of the occupation and the currents to the cotunneling process is different before and after the voltage pulse. Before the voltage pulse, the system is dominated by the cotunneling process because the Fermi surface of the both leads are nearly equivalent but is lower than the dot energy level (µ L ≃ µ R ≃ 0 < ǫ = 2Γ). The aligned fermi surfaces double the peak amplitude of the displacement current which also gives the zero net current, while the occupation has a corresponding change regarding the change of the displacement current. With the voltage increasing, the Fermi surface of the left lead moves up over the dot energy level while that of the right lead moves down below the dot energy level such that the system gradually approaches the sequential tunneling regime. This drives the electron flowing from the left lead to the dot and then to the right lead. Such a process results in a sensitive response of all physical quantities, the electron occupation in the dot, the displacement and the net currents. Then with the voltage decaying to zero, the electron residing in the dot is favorable to cotunnel to the left lead which gives a negative current and finally reaches the steady-state. The transient electron dynamics with the non-linear response to this Gaussian pulse is in particular useful for quantum feedback control to the electron states in the dot through the transient current that we will study in the future. It is also favorable to study the quantum capacitance and inductances. 72,76 Note that the above transient dynamics starts with a zero initial occupation, i.e., ρ 00 (t 0 ) = 1 and ρ 11 (t 0 ) = 1 − ρ 00 (t 0 ) = 0 [N (t 0 ) = 0]. This implies that the last term of Eq. (30) which is often ignored in the nonequilibrium Green function technique 11,17 has no contribution to the transient current in the above numerical calculations. If the dot is initially occupied, the transient current will be quite different although the steady-state limit is the same because the initial electron distribution vanishes in the steady-state limit. In Fig. 9, we show the exact numerical result with such a situation where we take the simple step-pulse voltage as an example. It shows that the initial occupation in the central region has measurable contributions in studying the transient dynamics, especially for the ultrafast (extremely short time) operations in a quantum device.
Combining all these analysis together, as one can see the exact numerical solutions presented here have demon- strated that both the electron states in the dot and the transient currents past through it have a clear non-linear response to the external bias voltage pulses, in particular within a short time scale after the pulse is turned on. These ultrafast non-linear response properties will provide very useful information for the manipulation of the device states as well as the quantum feedback controls for practical applications. However, since we only consider here a very simple device with a single-level dot, the quantum coherence and decoherence dynamics to the transient current do not be manifested in these numerical solutions. To demonstrate the quantum coherence properties in the transient transport dynamics, it is necessary to have the device containing at least two levels (or a single level with electron spin degrees of freedom) in the central region. Also, in the above numerical calculations, we take a rather large finite line width, W α = 20Γ. The non-Markovian memory structure will become more significant when the line width becomes narrower. 17, 21 We will examine in detail these features within the present theory in our future work.

V. SUMMARY AND PROSPECTIVE
In summary, we have established a nonequilibrium theory for the transient quantum transport dynamics via the Feynman-Vernon influence functional approach. We obtain an analytical relation between the master equation, Eq. (8), and the transient current, Eq. (19), determined by the same time-dependent (non-Markovian) coefficients. In particular, the master equation and the transient current are explicitly related to each other in terms of the superoperators acting on the reduced density matrix of the central system, see Eq. (23). The back-reaction effect of the gating electrodes to the central system is fully taken into account by the time-dependent coefficients in the master equation (8) through the dissipation-fluctuation integrodifferential equation (7). The non-Markovian memory structure is nonperturbatively built into the integral kernels in these equations of motion. The resulting transient current, Eq. (19) or (30), is rather simple and it recovers the steady-state current in the nonequilibirum Green function technique and the Landauer-Büttiker formula at the long time limit. This exact nonequilibrium formalism should provide a very intuitive picture of how the change of the electron quantum coherence in the central system is intimately related with the electron tunneling processes through the leads, and therefore non-linearly responses to the corresponding external bias controls. This theory is applicable to study a variety of quantum transport phenomena involving explicitly non-Markovian quantum decoherence behaviors, in both stationary and transient scenarios, at arbitrary temperatures of the different contacts in the weak Coulomb interaction regime.
As a simple illustration, we apply the theory to a simple model of electron tunneling though a single level quantum dot. We obtain all analytical solutions in the wide band limit (WBL) where the dependence of the initial electron occupation in the dot is explicitly manifested. Taking a more realistic spectral density with Lorentzian shape, we show that the Markov limit is a good approximation in the WBL. The non-Markovian memory effect is dominated by a finite line width of the spectral density. Under the ac bias voltage pulses (including the step pulse, the Gaussian pulse and the oscillation pulse), we have demonstrated the ultrafast nonlinear response of the electron occupation and currents to the ac bias. We find that the current evolutions are more sensitive to the energetic configuration of the dot than the occupation evolution. This feature is very useful for quantum feedback controls in quantum information processing. More applications will be presented in the future works. These include the decoherence transport dynamics in quantum dot devices such as quantum dot Aharonov-Bohm inteferometers; the nonequilibrium dynamics and the real time monitoring of spin polarization processes in nanostrucutres; also the transient transport dynamics in molecular electronics, as well as the application to bio-electronics such as DNA junctions, etc.
In the end, we should also point out that the present theory is developed without considering the electronelectron interaction in the central region and therefore it is mainly valid in the weak Coulomb interaction regime. It is not difficult to extend the present theory to the strong Coulomb Blockage regime by properly excluding the doubly occupied states in the central region of the nanostructure, as we have shown explicitly in [21]. Although studying the above two extreme limits, the extremely weak and the extremely strong Coulomb interaction regimes together, could reveal a significant understanding to various quantum transport phenomena, there is increasing discussion for the intermediate Coulomb interaction regime 77,78 where the analysis of transport dy-namics should become much more complicated. In this situation, the path integrals of Eqs. (3) and (14) may not be carries out exactly, and therefore it is not easy to find an exact master equation and an exact expression of the transient current. But the master equation, Eq. (8), and the transient current, Eq. (19), can still serve as a good approximation with respect to the saddle point approximation or loop expansion, 79 where the Coulomb interaction must be included self-consistently in the dissipation-fluctuation integrodifferential equation (7). This approximate treatment could provide a basic understanding to the quantum transport phenomena in the intermediate Coulomb interaction regime. The detailed extension of the present theory to the interacting electron systems is in progress now. Nevertheless, the analytical nonequilibrium theory we presented in this paper has the advantage of combining the decoherence dynamics with transient transport dynamics together to explore time-dependent physical phenomena and may also provide a basic theory for quantum feedback controls in various nanoelectronics devices.  As we see both the master equation (8) and the transient current (30) are completely determined by the propagating matrices of the stationary paths: u(τ ),ū(τ ) and v β (τ ). Here we shall show that these propagating matrices are directly related to the retarded, advanced and lesser Green functions in the nonequilibrium Green function technique. In our previous work, 21 the propagating matrices u(τ ),ū(τ ) and v β (τ ) were introduced to simplify the stationary path equations of motion (with the convention x = yz for x i = j y ij z j ) as follows: ξ(τ ) = u(τ )ξ(t 0 ) + v β (τ )[ξ(t) + ξ ′ (t)], (A.1a) ξ(τ ) + ξ ′ (τ ) =ū(τ )[ξ(t) + ξ ′ (t)]. (A.1b) In fact, Eq. (A.1) shows that u(τ ) is a propagating matrix of the forward stationary paths ξ(τ ) starting at t 0 , while v β (τ ) mixes the forward path ξ(τ ) and the backward path ξ ′ (τ ) started backwardly from t, andū(τ ) is a backward propagating matrix of the stationary paths. In fact, Eq. (A.1) shows that the transformation matrices u(τ ),ū(τ ) and v β (τ ) should be defined more precisely as u(τ ) ≡ u(τ, t 0 ),ū(τ ) ≡ u † (t, τ ) and v β (τ ) ≡ v β (τ, t): ξ(τ ) = u(τ, t 0 )ξ(t 0 ) + v β (τ, t)[ξ(t) + ξ ′ (t)], (A.2a) where the Grassmannian variables ξ(τ ) and ξ ′ (τ ) represent the forward and backward electron paths in the functional path integrals. Then Eq. (A.1) directly tells that u(τ, t 0 ) describes the electron propagation (represented by ξ(τ ) in the Grassmannian space) from the initial time t 0 to the time τ so that it is just the retarded Green function, namely, Eq. (33) is a justification for this relation. While, u † (t, τ ) describes the inverse propagation of the electron [or the backward propagation represented by ξ ′ (τ )] from the time t to the time τ such that it is indeed the advanced Green function: Meantime, Eq. (7a) indicates that the time correlation function of the α-reservoir: as a back-reaction effect of the α-lead to the central system, is the retarded self-energy. These relations are also justified by Eqs. (7a-7b). The function v β (τ, t) describes the electron propagation mixing the forward and backward paths so that it is related to the lesser Green function defined by G < ij (τ, t) ≡ i a † j (t)a i (τ ) in the nonequilibrium Green function formalism. In fact, it is not difficult to find the explicit solution of Eq. (7c):  where Σ < (τ 1 , τ 2 ) = ig β (τ 1 , τ 2 ), (A.8) is the lesser component of the self-energy. On the other hand, the single-particle reduced density matrix is related to the lesser Green function by ρ (1) (t) = −iG < (τ, t)| τ =t . From the relation of Eq. (29), we find that v β (τ ) can be expressed in terms of the lesser Green function as follows: =G r (τ, t 0 )G < (t 0 , t 0 )G a (t 0 , t) + τ t0 dτ 1 t t0 dτ 2 G r (τ, τ 1 )Σ < (τ 1 , τ 2 )G a (τ 2 , t). (A.9) This provides indeed a general solution of the lesser Green function in the nonequilibrium Green function technique. In the previous investigation of transient dynamics, the first term in Eq. (A.9) that sensitively depends on the initial electron distribution in the central regions are often ignored. 11,17,24 In fact, the second term in Eq. (A.9) can only be identified as the lesser Green function G < (τ, t) in the steady-state limit.
Using the above explicit relations between u(t),ū(τ ), v β (t) and the nonequilibrium retarded, advanced and lesser Green functions, we immediately obtain the transient current of Eq. (30) in terms of the nonequilibrium Green functions: dτ Tr Σ r α (t, τ )G < (τ, t) This current has exactly the same form obtained form the nonequilibrium Green function technique. 11,24 However, as we have pointed out in the practical applications, one usually uses the steady-state lesser Green function, namely ignoring the first term in Eq. (A.9). This term vanishes in the steady-state limit so that it will not affect the steady-state current. It can also be dropped if one takes the initial time t 0 → −∞ so that the central region is assumed to be in an empty state initially. However, this term which explicitly depends on the initial single particle reduced density matrix (including the initial electron occupation in each level and the electron quantum coherence between different levels in the central region) is crucial for practical manipulation of a real quantum device. Only in the wide band limit (WBL) where the non-local time correlation function g α (t, τ ) = iΣ r α (t, τ ) = Γ α δ(t − τ ), (A.11) the integral of the first term in Eq. (A.10) is reduced to the single particle reduced density matrix: G < (t, t) = iρ (1) (t)), which results in dτ Im e −i[ω(t−τ )+e R t τ dτ ′ Vα(τ ′ )] G a (τ, t) . (A.12) Thus, the ignored initial occupation dependence is recovered in the WBL.