Time-dependent quantum transport: A practical scheme using density functional theory

S. Kurth, G. Stefanucci, C.-O. Almbladh, A. Rubio, and E.K.U. Gross Institut für Theoretische Physik, Freie Universität Berlin, Arnimallee 14, D-14195 Berlin, Germany Solid State Theory, Institute of Physics, Lund University, Sölvegatan 14 A, S-22362 Lund, Sweden Departamento de Física de Materiales, Facultad de Ciencias Químicas, UPV/EHU, Unidad de Materiales Centro Mixto CSIC-UPV/EHU and Donostia International Physics Center (DIPC), San Sebastián, Spain (Dated: October 6, 2018)


I. INTRODUCTION
During the last decades, the size of electronic circuits has continuously been reduced. Today, systems like quantum wires and quantum dots are routinely produced on the nanometer scale. Recently, the seemingly ultimate limit of minituarization has been achieved by several experimental groups who were able to place single molecules between two macroscopic electrodes. 1,2 To describe transport properties on such a small scale, a quantum theory of transport is required. 3,4 A cornerstone of such a theory is the Landauer-Büttiker formalism 5,6 which provides a method to compute the steady-state current of non-interacting electrons for meso-or nanoscopic systems connecting two (or more) macroscopic electrodes.
Alternatively, the technique of nonequilibrium Green functions (NEG) 7,8 has been used to tackle quantum transport. Studies using the NEG approach typically use tight-binding-like model Hamiltonians to describe the combined system electrodes plus "device". A well known scheme is the one introduced by Caroli et al. 9,10 In the remote past the left and right electrodes are disconnected and in equilibrium at two different chemical potentials; the conducting part of the Hamiltonian is switched on adiabatically and eventually a steady-state develops. Within this contacting approach also time-dependent transport phenomena have been investigated. 11 Caroli et al. discussed non-interacting systems only. Their ap-proach has later been extended to account for short-range electron-electron interaction and for interaction with vibrations in the device region. 12 An excellent overview of the field can be found in the book by Haug and Jauho 13 and in Ref. 4. Despite its appeal, the above scheme has limitations since the time-dependent perturbation is the tunneling Hamiltonian, connecting the electrodes to the device, rather than the external electric field.
Cini proposed another scheme 14 also based on NEG. Here, the system electrodes plus "device" is connected and in equilibrium in the remote past. The timedependent perturbation is the external scalar potential. It has been shown 15 that for non-interacting systems the contacting approach and the Cini approach yield the same current in the long-time limit and that in the dc case the steady state current does not depend on the history of the applied potential. Moreover, the Cini scheme is well suited for a density functional extension since the electrons are driven out of equilibrium by a local potential rather than by a non-local one (see below).
With recent experimental progress to place single molecules as devices between macroscopic electrodes there also has been considerable activity to describe transport through these systems on an ab initio level. Most approaches are based on a self-consistency procedure first proposed by Lang. 16 In this steady-state approach based on density functional theory (DFT), exchange and correlation is approximated by the static Kohn-Sham (KS) potential and the charge density is obtained self-consistently in the presence of the steady current. However, the original justification involved subtle points such as different Fermi levels deep inside the left and right electrodes and the implicit reference of non-local perturbations such as tunneling Hamiltonians within a DFT framework. (For a detailed discussion we refer to Ref. 15.) The steady-state DFT approach has been further developed 17,18,19,20 and the results have been most useful for understanding the qualitative behavior of measured current-voltage characteristics. Quantitatively, however, the theoretical I-V curves often differ from the experimental ones by several orders of magnitude. 21 Several explanations are possible for such a mismatch: models are not sufficiently refined, parasitic effects in measurements have been underestimated, the characteristics of the molecule-contact interface are not well understood and difficult to address given their atomistic complexity. Adding to the theoretical reason for this discrepancy is the fact that the transmission functions computed from static DFT have resonances at the non-interacting Kohn-Sham excitation energies which in general do not coincide with the true excitation energies. Furthermore, different exchange-correlation functionals lead to DFT-currents that vary by more than an order of magnitude. 22 Excitation energies of interacting systems are accessible via time-dependent (TD) DFT. 23,24 In this theory, the time-dependent density of an interacting system moving in an external, time-dependent local potential can be calculated via a fictitious system of non-interacting electrons moving in a local, effective time-dependent potential. Therefore this theory is in principle well suited for the treatment of nonequilibrium transport problems. 25 A basic issue is that most exchange-correlation functionals have been derived under equilibrium conditions and their application to non-equilibrium problems should be analyzed in more detail. However, this is beyond the scope of the present work.
Before a TDDFT calculation of transport can be tackled, a number of technical problems have to be addressed. In particular, one needs a practical scheme for the propagation of the time-dependent Schrödinger equation for an infinitely large system. Of course, since one can in practice only deal with finite systems this can only be achieved by applying the correct boundary conditions. The problem of so-called "transparent boundary conditions" for the time-dependent Schrödinger equation has been attacked by many authors. For a recent overview, the reader is referred to Ref. 26.
In this paper we present a propagation scheme which is particularly designed to be used for the calculation of time-dependent transport problems. In Section II, we combine the Cini scheme with TDDFT and we develop a general formalism based on the propagation of Kohn-Sham orbitals in open systems. In Section III we will address the question of how to obtain the correct initial states for the propagation. An algorithm for the timeevolution of open systems is proposed in Section IV. It is based on a modified version of the Crank-Nicholson algorithm. Section V describes some details of our numerical implementation and Section VI gives results for several one-dimensional model systems. We draw our conclusions in Section VII.

II. GENERAL FORMULATION
We consider an electrode-junction-electrode system which is initially in equilibrium (t < 0). The system is contacted and no current flows through the junction, the charge density of the electrodes being perfectly balanced, see Fig. 1. Therefore, the system initially is in its ground state which, due to the Hohenberg-Kohn theorem, 27 is a functional of the density. This density can then be computed in the usual way by n(r, 0) = occ |ψ s (r, 0)| 2 where the sum is over the occupied Kohn-Sham orbitals ψ s (r, 0), i.e., the eigenfunctions of the Kohn-Sham Hamiltonian H(0) with eigenergy below the Fermi energy. Here and in the following we use boldface notation to denote operators in one-electron Hilbert space.
To observe a current we drive the system out of equilibrium by exposing the electrons to an external electric potential (bias). Without loss of generality we will assume that this external bias vanishes for times t ≤ 0. According to the Runge-Gross theorem 23 , the time-dependent density can be calculated by evolving the KS orbitals according to the KS equation of TDDFT iψ s (t) = H(t)ψ s (t) where H(t) is the time-dependent KS Hamiltonian. Thus, n(r, t) = occ |ψ s (r, t)| 2 and the continuity equation is the KS current density. Due to the Runge-Gross theorem and the continuity equation one can deduce that the longitudinal part of the KS current density equals the longitudinal part of the true current density. This need not be true for the transverse part of the current density. However, the transverse part of the current density does not contribute to the total current which can then be calculated by a surface integral wheren is the unit vector perpendicular to the surface element dσ and the surface S is perpendicular to the longitudinal geometry of our system. In order to propagate the KS orbitals we need to solve the Schrödinger equation for a macroscopic and non-periodic system. This goal is hopeless unless we know the dynamics of the remote parts of the system. We restrict ourselves to metallic electrodes. Then, the external potential and the disturbance introduced by the device region are screened deep inside the electrodes. As the system size increases, the remote parts are less disturbed by the junction and the density inside the electrodes approaches the equilibrium bulk-density. Thus, the macroscopic size of the electrodes leads to an enormous simplification since the initial-state self-consistency is not disturbed far away from the constriction.
It is convenient to partition the system into three main regions: a central region C consisting of the junction and a few atomic layers of the left and right electrodes and two regions L, R which describe the left and right bulk electrodes. Only the central device region C will be treated explicitly. Our scheme accounts for the full dynamical screening in the central region. It can be further refined by taking into account screening effects also deeper in the electrodes at the level of linear response, with a limited increase in numerical efforts. (These effects might be of importance in the initial transient phase where long-range plasma oscillations in the electrodes may occur.) According to the above partitioning, the original KS Hamiltonian can be written as a 3 × 3 block matrix, and the Schrödinger equation reads where ψ α is the projected wave-function onto the region α = L, R, C. We can solve the diffential equation for ψ L and ψ R by introducing the retarded Green function g α for electrode α, which satisfies i d dt − H αα (t) g α (t, t ′ ) = δ(t − t ′ ), α = L, R, with boundary conditions g α (t + , t) = −i and g α (t, t + ) = 0. Then, we have for α = L, R Using Eq. (3), the equation for ψ C can be written as where Σ = α=L,R H Cα g α H αC is the self-energy which accounts for the hopping in and out of region C. Thus, for any given KS orbital we can evolve its projection onto the central region by solving Eq. (4) in region C. Eq. (4) has also been derived elsewhere (for static Hamiltonians). 28 To summarize, all the complexity of the infinite electrode-junction-electrode system shown in Fig. 1 has been reduced to the solution of an open quantum-mechanical system (the central region C) with proper time-dependent boundary conditions.

III. COMPUTATION OF EXTENDED EIGENSTATES
Eq. (4) of the previous Section is the central equation of our approach to time-dependent transport. It is a reformulation of the original time-dependent Schrödinger equation (2) of the full system in terms of an equation for the central (device) region only. The coupling to the leads is taken into account by the lead Green functions g α , α = L, R. Eq. (4) has the structure of a timedependent Schrödinger equation with two extra terms: the term involving the self energy Σ which we will also denote as the memory integral since it involves the wavefunction in the central region at previous times during the propagation. 29 The second term describes the injection of particles induced by a non-vanishing projection of the initial wave-function onto the leads.
Eq. (4) is first order in time, therefore we need to specify an initial state which is to be propagated. We want to study the time evolution of systems perturbed out of their equilibrium ground state. Of course, the ground state of our noninteracting system is the Slater determinant of the occupied eigenstates of the full, extended Hamiltonian in equilibrium, H(t < 0) = H s . The practical question then arises how one can obtain these eigenstates without having to deal explicitly with the extended Hamiltonian. Note that, unlike in a bulk solid, the translational symmetry is broken in the present device situation.
In the present Section we propose a solution of this problem which is based on the partitioning approach used in many steady-state transport calculations (see, e.g., Ref. 31). The retarded Green function of the static Hamiltonian in the energy domain is determined by The Green function G(E) of the full system can be written in the same block structure as the Hamiltonian Eq. (5) can be solved for the block of the Green function referring only to the central region with the retarded Green function of lead α (8) and the unit matrix 1 α in region α. This Green function enters as a central ingredient into the Fisher-Lee relation 32 for the calculation of the transmission function. Through the coupling to the leads it provides for level broadening of the isolated central part, but it also contains information on the eigenstates of the extended system.
In order to illustrate the central idea of our method to extract the extended eigenstates from the Green function we consider H s to be the discretized form of a continuous HamiltonianĤ s (r). The continous Green function and the discretized one for uniform lattice spacing ∆x are connected by where N = 1, 2, 3 is the number of spatial dimensions of the problem. We choose the convention that a singleparticle orbital ψ Ej of the HamiltonianĤ s is uniquely specified by its eigenenergy E and a label j for the d E degenerate orbitals of this energy. Using the Lehmann representation and assuming thatĤ s is invariant under time-reversal, the imaginary part of G is (10) Multiplying Eq. (10) by ψ El (r ′ ), integrating over r ′ , using the orthogonality of the single particle states and dividing by the density of states Eq. (11) has the structure of an eigenvalue equation where the energy eigenstate ψ El is also an eigenstate of the integral operator −Im [G(r, r ′ , E)]/(πD(E)) with eigenvalue γ(E). For a given energy E, this integral operator has d E different degenerate eigenstates. We note that Eq. (10) is valid for all points in space, in particular also for both r and r ′ representing points in the central region. In this case we know the Green function G CC through Eqs. (7) and (9). Below we show that the eigenfunctions of Im[G CC ] can be expressed as linear combination of the ψ El . Let us consider the matrix formed by the elements where the integration is over the central region only. This matrix is Hermitian and can be diagonalized, i.e., with λ j real. Next we compute the matrix elements of the Green function with respect to the functions After a straightforward manipulation one finds which shows explicitly that the functions a Ej (r) diagonalize Im [G CC ] in the central region and that the eigenvalues are positive. Since any linear combination of degenerate eigenstates is again an eigenstate, diagonalizing Im [G CC (E)] gives us one set of linearly independent, degenerate eigenstates of energy E. In our practical implementation described in more detail in Section V, we diagonalize where is the total density of states in the central region. If we use N g grid points to describe the central region, the diagonalization in principle gives N g eigenvectors but only a few have the physical meaning of extended eigenstates at this energy. It is, however, very easy to identify the physical states by looking at the eigenvalues: only few eigenvalues (for the simple examples we studied either one or two) are nonvanishing and they always add up to unity. The corresponding states are the physical ones. All the other eigenvalues are zero (or numerically close to zero) and the corresponding states have no physical meaning. The procedure described above gives the correct extended eigenstates only up to a normalization factor. When diagonalizing Eq. (17) with typical library routines one obtains eigenvectors which are normalized to the central region. Physically this might be incorrect. Therefore, the normalization has to be fixed separately. In the example of Section V we fixed the norm by matching the wavefunction for the central region to the known form (and normalization) of the wavefunction in the macroscopic leads.
It should be emphasized that the procedure described here for the extraction of eigenstates of the extended system from G CC (E) only works in practice if E is in the continuous part of the spectrum due to the sharp peak of the delta function in the discrete part of the spectrum.
Eigenstates in the discrete part of the spectrum can be found by considering the original Schrödinger equation for the full system: Using again the block structure of the Hamiltonian this can be transformed into an effective Schrödinger equation for an energy-dependent Hamiltonian for the central region only: This equation has solutions only for certain values of E which are the discrete eigenenergies of the full Hamilto-nianĤ s . Therefore, one can find these states by iteration. We have succesfully tested this idea for systems where the analytic solutions are known. However, since the main focus of the present work is transport where the continuum states are the essential physical ingredient, we will not deal with the states in the discrete spectrum for the remainder of this paper. Those states might play a role in the description of charge-accumulation in molecular transport, as, e.g., in Coulomb blockade phenomena.

IV. ALGORITHM FOR TIME EVOLUTION
In order to calculate the longitudinal current in an electrode-junction-electrode system we need to propagate the Kohn-Sham orbitals. The main difficulty stems from the macroscopic size of the electrodes whose remote parts,ultimately, are taken infinitely far away from the central, explicitly treated, scattering region C.
The problem can be solved by imposing transparent boundary conditions 26 at the electrode-junction interfaces. Efficient algorithms have been proposed for wavepackets initially localized in the scattering region and for Hamiltonians constant in time. In this Section we propose an algorithm well suited for delocalized initial states, as well as for localized ones, evolving with a timedependent Hamiltonian.
Let H(t) be the time-dependent KS Hamiltonian. We partition H(t) as in Section II. The explicitly treated region C includes the first few atomic layers of the left and right electrodes. The boundaries of this region are chosen in such a way that the density outside C is accurately described by an equilibrium bulk density. It is convenient to write H αα (t), with α = L, R, as the sum of a term H s αα which is constant in time and another term U α (t) which is explicitly time-dependent, H αα (t) = H s αα + U α (t). In configuration space U α (t) is diagonal at any time t since the KS potential is local in space. Furthermore, the diagonal elements U α (r, t) are spatially constant for metallic electrodes. Thus, and In this way, the only term inH(t) that depends on t is H CC (t). For any given initial state ψ(0) = ψ (0) we calculate ψ(t m = m∆t) = ψ (m) by using a generalized form of the Cayley method ] and δ = ∆t/2. It should be noted that our propagator is norm conserving (unitary) and accurate to second-order in δ, as is the Cayley propagator. 33 Denoting by ψ α the projected wave function onto the region α = R, L, C, we find from Eq. (23) Here, H ]. The source term S (m) describes the injection of density into the region C, while the memory term M (m) is responsible for the hopping in and out of the region C. In terms of the propagator for the uncontacted and undisturbed α electrode the source term can be written as with For a wave packet initially localized in C the projection onto the left and right electrode ψ (0) α vanishes and S (m) = 0 for any m, as it should be. The memory term is more complicated and reads where The quantities Q (m) α depend on the geometry of the system and are independent of the initial state ψ (0) .
Below we propose a recursive scheme to calculate the Q (m) α 's for those system geometries having semiperiodic electrodes along the longitudinal direction, see Fig. 2. In this case H s αα has a tridiagonal block form where h α describes a convenient cell and V α is the hopping Hamiltonian between two nearest neighbor cells.
Without loss of generality we assume that both h α and V α are square matrices of dimension N α × N α . Taking into account that the central region contains the first few cells of the left and right electrodes, the matrix Q (m) α has the following structure

FIG. 2:
Schematic sketch of an electrode-junction-electrode system with semiperiodic electrodes. and are given by where the subscript (1, 1) denotes the first diagonal block of the matrix in the square brackets. We introduce the generating matrix function which can also be expressed in terms of continued matrix fractions (see Appendix ) The q (m) α 's can be obtained from From Eqs. (35) and (34) one can build up a recursive scheme. Let us define p −1 α (x, y) = x + iyδh α + y 2 δ 2 q α (x, y) and p Using the identity 1 with p This concludes the description of our algorithm for the propagation of the time-dependent Schrödinger equation for extended systems. It is worth mentioning an additional complication here which arises for the propagation of a time-dependent Kohn-Sham equation. This complication stems from the fact that in order to compute ψ (m+1) C at time step m + 1 one needs to know the timedependent KS potential at the same time step which, via the Hartree and exchange-correlation potentials, depends on the yet unknown orbitals ψ (m+1) C . Of course, the solution is to use a predictor-corrector approach: in the first step one approximates H

V. IMPLEMENTATION DETAILS
All the methodological discussion above is general and can be applied to general device configurations as long as they can be mapped into a longitudinal-like geometry as described in Fig. 2. In order to demonstrate the feasibility of the scheme described in the previous Sections we have implemented it for one-dimensional model systems. The extension to real molecular-device configurations is presently under development. 34 We have used a simple three-point discretization for the second derivative with equidistant grid points x i , i = 1, . . . , N g and spacing ∆x. Within this approximation matrices of the form H Cα M H αC which are N g × N g matrices and appear, e.g., in Eq. (7) or (29), have only one nonvanishing matrix element. For α = L this is the (1, 1) element, for α = R it is the (N g , N g ) element.
In order to proceed we have to specify the nature of the leads and therefore the lead Green function. Here we choose the simplest case of semi-infinite, uniform leads at constant potential U α0 . In this case, the Green function (8) in the energy domain can be given in closed form: The abscissa x α0 is the position of the interface between the lead and the device region and x k = x α0 ± k∆x, where the plus sign applies for α = R and the minus sign for α = L.
The results of the procedure for calculating extended eigenstates as described in Section III is illustrated in Fig. 3 for a square potential barrier with zero potential in both leads. In the upper panel we have the square modulus of eigenstates at an energy below the barrier height while in the lower panel eigenstates with energy higher than the barrier are shown. The states result from diagonalization of Eq. (17). In order to obtain the normalization constant we compute the logarithmic derivative at the boundary of the central region numerically and match it to the analytic form in the lead to obtain the phase shift δ α : where q = 2Ẽ α . Knowing the phase shift we can rescale the wavefunction such that it matches with the analytic form sin(q(x − x α0 ) + δ α ) at the interface. Of course, this form of the extended states only applies forẼ α > 0 but as long as E is in the continuous part of the spectrum, it is correct at least for one of the leads. This is sufficient to determine the normalization. The states obtained numerically with this procedure coincide with the known analytical results. We then implemented the propagation scheme presented in the previous Section. Within our three-point Although the quadratic equation has two solutions, the above choice for q (0) α is dictated by the fact that the Taylor expansions for small δ of Eqs. (41) and (34) have to coincide. Using this result we then solved the iterative scheme to obtain the q (m) α for m ≥ 1. As a first check on the propagation method we propagated a Gaussian wavepacket which, at initial time t = 0, is completely localized in the central device region. (The source terms S (m) then vanish identically). As can be seen in Fig. 4, the wavepacket correctly propagates through the boundaries without any spurious reflections.
For the propagation of the extended initial states (the eigenstates of the unperturbed system) we also need to implement the source terms S (m) . In the following we assume that the left and right leads are at the same potential initially so that the analytic form of the initial states is in both leads given by sin(q(x − x α0 ) + δ α ) = cos(δ α ) sin(q(x − x α0 )) + sin(δ α ) cos(q(x − x α0 )). The propagation of the term proportional to sin(q(x − x α0 )) is trivial since this is an eigenstate of the lead Hamiltonian with energy ε q = q 2 /2. Therefore, if -in discretized form -ψ and similarly for the left lead. Here, H s RR is the static part of the right-lead Hamiltonian, g R the corresponding Green function according to definition (25) and e R = (0, . . . , 0, 1) T is a unit vector. The propagation of the part proportional to cos(q(x − x α0 )) is more complicated since this is not an eigenstate of the lead Hamiltonian. We define the function R(x, y) from where ψ and using one arrives at where q R (x, y) is given by Eq. (33). Defining one finds and finally One may proceeds along the same lines for extracting the left component.
To test our implementation we have propagated eigenstates of the extended system. As expected, these states just pick up an exponential phase factor exp(−iEt) during the propagation.

VI. EXAMPLES
We are now in a position to apply our algorithm to the calculation of time-dependent currents in onedimensional model systems. The systems are initially in thermodynamic equilibrium. At time t = 0, a bias is switched on in the electrodes.
As a first example we considered a system where the electrostatic potential vanishes identically both in the left and right leads as well as in the central region which is explicitly propagated. Initially, all single particle levels are occupied up to the Fermi energy ε F . At t = 0 a constant bias is switched on in the leads and the timeevolution of the system is calculated. We chose the bias in the right lead as the negative of the bias in the left lead, U R = −U L . The current is calculated from Eq. (1): where the prefactor 2 comes from spin and k F = √ 2ε F is the Fermi wavevector.
The numerical parameters are as follows: the Fermi energy is ε F = 0.3 a.u., the bias is U L = −U R = 0.05, 0.15, 0.25 a.u., the central region extends from x = −6 to x = +6 a.u. with equidistant grid points with spacing ∆x = 0.03 a.u.. The k-integral in Eq. (50) is discretized with 100 k-points which amounts to a propagation of 200 states. The time step for the propagation was ∆t = 10 −2 a.u.
In Fig. 5 we have plotted the current densities at x = 0 as a function of time for different values of the applied bias. As a first feature we notice that a steady state is achieved, in agreement with the results of Ref. 15. The corresponding steady-state current I can be calculated from the Landauer formula. For the present geometry this leads to the steady current . We see that our algorithm yields the same answers. Second, we notice that the onset of the current is delayed in relation to the switching time t = 0. This is easily explained by the fact that the perturbation at t = 0 happens in the leads only, e.g., for |x| > 6 a.u., while we plot the current at x = 0. In other words, we see the delay time needed for the perturbation to propagate from the leads to the center of our device region. We also note that the higher the bias the more the current overshoots its steady-state value for small times after switching on the bias. Finally it is worth mentioning that increasing the bias not necessarily leads to a larger steady-state current.
In the second example we studied a double square potential barrier with electrostatic potential V (x) = 0.5 a.u. for 5 a.u. ≤ |x| ≤ 6 a.u. and zero otherwise. This time we switch on a constant bias in the left lead only, i.e., U R = 0. The Fermi energy for the initial state is ε F = 0.3 a.u.. The central region extends from x = −6 to x = +6 a.u. with a lattice spacing of ∆x = 0.03 a.u.. Again, we use 100 different k-values to compute the current and a time step of ∆t = 10 −2 a.u..
In Fig. 6 we plot the current at x = 0 as a function of time for several values of the applied bias U = U L . Again, a steady state is achieved for all values of U . As discussed in Fig. 5 the transient current can exceed the steady current; the higher the applied voltage the larger is this excess current and the shorter is the time when it reaches its maximum. Furthermore, the oscillatory evolution towards the steady current solution depends on the bias. For higher bias the frequency of the transient oscillations increases. For small bias the electrons at the bottom of the band are not disturbed and the transient process is exponentially short. On the other hand, for a bias close to the Fermi energy the transient process decays as a power law, due to the band edge singularity. As pointed out in Ref. 15, for non-interacting electrons the steady-state current develops by means of a pure dephasing mechanism. In our examples the transient process occurs in a femtosecond time-scale, which is much shorter than the relaxation time due to electron-phonon interactions.
In Fig. 7 we plot the time evolution of the total number of electrons in the device region for the same values of U L . We see that as a result of the bias a quite substantial amount of charge is added to the device region. This result has important implications when simulating the transport through an interacting system as the effective (dynamical) electronic screening is modified due not only to the external field but also to the accumulation of charge state in the molecular device. This illustrates that linear response might not be an appropriate tool to tackle the dynamical response and that we will need to resort to a full time-propagation approach as the one of the present work. Here we emphasize that all our calculations are done without taking into account the electron-electron interaction. If we had done a similar calculation with the interaction incorporated in a time-dependent Hartree or time-dependent DFT framework we would expect the amount of excess charge to be reduced significantly as compared to Fig. 7.
In Fig. 8 we show time-dependent currents for the same double barrier as in Fig. 6 for two different ways of applying the bias in the left lead: in one case the constant bias U 0 is switched on suddenly at t = 0 (as in Fig. 6), in the other case the constant U 0 is achieved with a smooth switching U (t) = U 0 sin 2 (ωt) for 0 < t < π/(2ω). As expected and in agreement with the results of Ref. 15, the same steady state is achieved after the initial transient time. However, the transient current clearly depends on how the bias is switched on.
In the final example we address the simulation of actransport. We computed the current for a single square potential barrier with V (x) = 0.6 for |x| < 6 and zero otherwise. Here we applied a time-dependent bias of the form U L (t) = U 0 sin(ωt) to the left lead. The right lead remains on zero bias. The numerical parameters are: Fermi energy ε F = 0.5 a.u., device region from x = −6 to x = +6 a.u. with lattice spacing ∆x = 0.03 a.u.. The number of k-points is 100 and the time step is ∆t = 10 −2 a.u.. In Fig. 9 we plot the current at x = 0 as a function of time for different values of the parameter U 0 = 0.1, 0.2, 0.3 a.u. The frequency was chosen as ω = 1.0 a.u. in both cases. Again, as for the dc-calculation discussed above, we get a transient that overshoots the average current flowing through the constriction; again, this excess current is larger the higher the applied voltage. Also, after the transient we obtain a current through the system with the same period as the applied bias. Note, however, that (especially for the large bias), the current is not a simple harmonic as the applied ac bias.

VII. SUMMARY AND OUTLOOK
In the present work we have presented a formally rigorous approach towards the description of charge transport using an open-boundary scheme within TDDFT. We have implemented a specific time-propagation scheme that incorporates transparent boundaries at the device/lead interface in a natural way. In order to have a clear definition of a device region in Fig. 1 we assumed that an applied bias can be described by adding a spatially constant potential shift in the macroscopic part of the leads. This implies an effective "metallic screening" of the constriction. The screening limits the spatial extent of the induced density created by the bias potential or the external field to the central region. Our treatment can be further refined to include dynamical screening deep inside the electrodes on the level of linear response, which might be of importance for the initial transient. Our time-dependent scheme allows to extract both ac and dc I/V device characteristics and it is ideally suited to describe external field (photon) assisted

processes.
In order to illustrate the performance of the method, we have implemented it for one-dimensional models and we have shown: i) How to extract the proper initial extended states to be propagated. ii) How to incorporate perfect transparent boundaries for the time propagation. iii) A steady-state current is always reached upon application of a dc bias. The transient process occurs on a time-scale much shorter than the relaxation time due to electron-phonon interaction. In the case of systems without any source of dissipation it is known that the steady-current is independent of the history of the process. 15 We have explicitly demonstrated this history independence for two different switching processes of the external bias. However, if we allow for dissipation either through electron-electron or electron-phonon interactions, the current versus voltage characteristics might depend on the history. For instance, hysteresis loops due to different transient electronic/geometrical device configurations are possible. This effect will be more dramatic in the case of ac-driving fields of high frequencies where the system may not have enough time to respond to the perturbation. iv) A periodic ac current is reached upon perturbation with a monochromatic field.
Previous work on time-dependent quantum transport mainly uses the idea of Caroli. 9,10 This approach is at the core of the Landauer derivation and has the problem of using different chemical potentials in different parts of the system. This implies that the initial state is not a ground state of the entire, contacted system. Furthermore, the time-dependent perturbation is a tunneling Hamiltonian which is nonlocal in space. Therefore, it cannot be combined with TDDFT since there the time-dependent potential is local. Here, we have used the scheme of Cini 14 which can be combined with TDDFT: We start the calculation from a well-defined thermodynamic equilibrium configuration, therefore the scheme is thermodynamically consistent. Then we apply an external potential that in general is time dependent. By virtue of the Runge-Gross theorem, the time-evolution of this quantum system is completely determined by the knowledge of the time-dependent density. In the steady state regime the occupation of left and right moving carriers is dictated by the time-dependent Schrödinger equation.
Another thermodynamically consistent scheme has been put forward by Kamenev and Kohn. 35 They used the microscopic quantum-kinetic formulation of conductivity and worked with the Kubo formula 36 in the linear time-dependent Hartree regime. They consider a closed system (ring) where the current in the system results from an external driving vector potential. This approach also overcomes the problem of having two chemical potentials. They have shown that it is possible to recover the Landauer result, i.e., the universal quantum of conductance 2e 2 /h independent of the length and material. Since the Kamenev-Kohn approach uses a vector potential rather than a scalar potential, time-dependent current DFT rather than TDDFT would be the natural density-functional extension.
Most theoretical approaches to transport in molecular electronic devices adopt open boundary conditions and assume that transport is ballistic, i.e., under steady state conditions inelastic collisions are absent from a molecular structure and its contacts. 4,5,6,31 Dissipation occurs only in the idealized macroscopic reservoirs connected by leads to the molecular device. This central role of the reservoirs in the process of dissipation is a valid approximation, particularly when the applied bias is small and a device operates in a linear regime. When inelastic scattering dominates this picture is not applicable. In particular, experiments 37,38 indicate that inelastic scattering with lattice vibrations is present at sufficiently large bias, causing local heating of contacts and molecular devices. Energy transfer to the lattice may also cause atomic migration and result in dramatic changes in the device characteristics (also they may give phonon replica structure in the measured conductance). The modelling of a many-electron system out of equilibrium coupled to lattice vibrations is a real theoretical challenge. 39 Electron correlations are also important in molecular conductors, for example, Coulomb blockade effects dominate the transport in quantum dots. Short-range electron correlations seems to be relevant in order to get quantitative description of I/V characteristics in molecular constrictions. 40,41,42 In particular it is commonly assumed that the energy scales for electron-electron and electronphonon interactions are different and could be treated separately. However, the metallic screening of the electrodes considerably reduces the Coulomb-charging energy (from eV to meV scale). In this regime the energy scale for the two interactions merge and they need to be treated on the same footing posing some additional theoretical challenge.
Other approaches put forward in the literature directly look for a homogeneous current-carrying state either based on a a maximum entropy principle 43 or by a imposing the current through Lagrange-multipliers. 44 In those approaches it is implicitly assumed that the origin of the homogeneous current is independent whether it is introduced by reservoirs or by external fields. This is indeed the case for independent electrons but once dissipation is built in the system might exhibit a dependence on the history of the applied bias,(e.g., possible appearance of hysteresis loop in the current versus voltage characteristics).
It is clear that the quality of the TDDFT functionals is of crucial importance. In particular, exchange and correlation functionals for the non-equilibrium situation are required. Time-dependent linear response theory for dc-steady state has been implemented in Ref. 45 within TDLDA assuming jellium-like electrodes (mimicked by complex absorbing/emitting potentials). It has been shown that the dc-conductance changes considerably from the standard Landauer value. Therefore, a systematic study of the TDDFT functionals themselves is needed. A step beyond standard adiabatic approximations and exchange-only potentials is to resort to manybody schemes as recently done for the characterization of optical properties of semiconductors and insulators. 46 Another path is to explore in depth the fact that the true exchange-correlation potential is current dependent. 47,48 An appealing feature of the present approach is that electron-electron and electron-ion correlations and dissipation would in principle be described correctly in twocomponent TDDFT. 49 At present we are implementing our propagation scheme for real 3D systems 34 within the real-space realtime TDDFT code, octopus. 50 We are also exploring the possibility of a semiclassical description of the electronion coupling in order to avoid the complexity of multicomponent DFT and the problems related to mixed quantum classical approaches (i.e., Ehrenfest dynamics) as they fail to describe the long-term inelastic electronphonon scattering correctly. 39 where Q 00 is the first N × N block of M −1 0 . It is now convenient to introduce the matrix M n obtained by dropping the first nN lines and nN columns of M 0 . Then, in terms of the rectangular matrixB n = [B n , 0, 0, . . .], we have