Non-equilibrium current and electron pumping in nanostructures
We discuss a numerical method to study electron transport in mesoscopic devices out of equilibrium. The method is based on the solution of operator equations of motion, using efficient Chebyshev time propagation techniques. Its peculiar feature is the propagation of operators backwards in time. In this way the resource consumption scales linearly with the number of states used to represent the system. This allows us to calculate the current for non-interacting electrons in large one-, two- and three-dimensional lead-device configurations with time-dependent voltages or potentials. We discuss the technical aspects of the method and present results for an electron pump device and a disordered system, where we find transient behaviour that exists for a very long time and may be accessible to experiments.
Electron transport through mesoscopic devices contacted to leads is intensely studied in chemistry and physics (see e.g. Refs. [1, 2]). The conceptual basis for theoretical studies is the non-equilibrium Green function (Keldysh) formalism . The Meir-Wingreen-formula  allows for the calculation of steady state currents. For time-dependent potentials or gate voltages the calculation of a current is a serious problem already for non-interacting electrons. One crucial point in most approaches is that the resource consumption scales quadratically with the number of system sites, since each Green function or operator product has two site indices. The necessity to study large systems conflicts with the rapid growth of computational demands, especially if long leads with long recurrence time are required. Recent studies therefore addressed mainly one-dimensional (1D) situations, e.g. by time-propagation of operator expectation values or of several single-particle eigenstates [5, 6, 7].
In the present contribution we numerically solve the equation of motion for a single fermion operator. In contrast to previous studies, the initial time coordinate is propagated backwards in time. Then the computational effort, in particular the memory consumption, scales linearly with the number of lattice sites used to represent device and leads. This allows us to calculate the non-equilibrium current even for large systems that could not be treated otherwise.
We consider the following situation: A device of -sites is contacted to two long leads extending along the -direction (see figure 1). This slab geometry includes the 2D () and 1D () case. For the kinetic energy in the Hamiltonian
we assume nearest-neighbour hopping () along each axis , and set as the unit of energy. The potential energy terms describe the voltage bias applied between the left/right lead, and the local potentials in the device. The current along a bond in the -direction is given by the expectation value of the corresponding operator
To specify the initial conditions, we assume that for times the device is isolated and uncharged, and the leads are in equilibrium at zero temperature (this situation can be described by setting along the contacts at and letting in ). At time device and leads are brought into contact, and electrons can flow onto the device. We then ask for the current through a section of the slab.
Time propagation of the many-fermion Fock state is not feasible, and the reduction of this problem to the propagation of several single-fermion states is restricted to non-interacting electrons and introduces additional errors that must be controlled. It appears to be more natural to propagate the current operators themselves. With the time evolution operator , defined by , , operators can be transformed to the Heisenberg picture
and expectation values are obtained from and the initial state as
Therefore we must represent in terms of operators at time ,
which allows evaluation of equation (4) since all initial expectation values are known.
We observe that two ways exist to obtain . The standard way starts from
In this equation, the commutator is given as a combination of operators , which is evaluated at time in the Heisenberg picture. Therefore solution of equation (6) requires time propagation of all operators . The total number of coefficients in the corresponding expansions equation (5) for these operators grows quadratically with the number of sites.
To avoid quadratic growth we here proceed the opposite way, starting from . In contrast to the previous case, we do not need to transform the commutator in
to the Heisenberg picture. Instead, it is only the single operator whose expansion
must be propagated in time. We do not need to keep track of all the other operators . Therefore memory consumption and computational effort are proportional to the number of sites. Making use of , each fermion operator in equation (2) can be propagated separately to obtain .
The price to be paid here is that equation (7) leads to propagation backwards in time, evolving the second time coordinate from the initial to the final . The main reason why we nevertheless consider this procedure is that, avoiding quadratic growth of memory consumption, only this way allows for the treatment of large systems, which are inaccessible in the first way. With modern Chebyshev techniques , time propagation itself is fast and accurate, and the additional effort for backward propagation is compensated by the reduction of the problem size. Note that for periodic time dependence, operators can be propagated iteratively by a full period.
Turning to the slab geometry, the eigenstates of decoupled infinite leads are discrete along the -, -direction and continuous along the -direction. In the actual calculation, infinite leads are replaced by long leads of length . The recurrence time revealing the artificial discretization is of the order , and must be larger than the maximal time in the calculations. Using operators for lead eigenstates to energy and operators for device sites as the in equation (8) the initial conditions specify the expectation values at time ,
where is the Fermi energy. Note that only these conditions change for finite temperatures.
In our first example in figure 2, a travelling potential wave
pumps electrons in a 1D system with zero voltage bias from the left to right (i.e. along the positive -direction). In the beginning electrons flow onto the initially uncharged device. The initial transient behaviour evolves into a ‘pseudo’ steady state (left panel), which persists over periods. The true steady state, with equal current at both sites of the device, is approached only in the extreme long time limit (right panel). Comparison to the case of a static potential, with the same but , shows that the ‘pseudo’ steady state persists on much longer time scales than the initial transient behaviour. This may allow for the experimental observation of deviations from the steady state, which here occur over periods up to the maximal propagation time. To allow for such long propagation time in our calculation, we choose long leads with . The full system has sites.
In our second example in figure 3, a constant voltage bias transports electrons through a disordered 3D device, with (static) random potentials with uniform probability distribution . A steady state is approached only for the short system with large bias . Longer systems () support localized states, whose presence causes oscillations in . These become more pronounced for smaller bias . The net current through the device is close to zero, since localized states do not contribute to transport. Note that the average current for is much smaller for than for , while for the system it is of the same magnitude.
In the given examples, up to lattice sites are treated for device and leads. Because of linear scaling in our method we must keep only that many elements in memory, and all calculations can be performed on standard desktop computers. With quadratic scaling, we would have to store about elements, requiring already GByte for complex numbers in double precision.
In conclusion, the numerical method discussed here allows for the calculation of time-dependent currents in large systems of non-interacting electrons, which are inaccessible to most other existing techniques. The propagation of operators backwards in time is a tolerable disadvantage in such cases. Apart from these achievements, the development of advanced solution techniques for complicated non-equilibrium Green function equations of motion appears more promising for future progress, especially with respect to the inclusion of dissipation or interaction. How such developments are possible within the context of the Chebyshev Space method  and Sparse Polynomial Space approach  will be discussed elsewhere. Data obtained with the present method serve as a reference to validate these new techniques.
The authors acknowledge support by Deutsche Forschungsgemeinschaft through SFB 652.
-  Cuniberti G, Fagas G and Richter K, eds 2005 Introducing molecular electronics (Springer)
-  Datta S 1995 Electronic Transport in Mesoscopic Systems (Cambridge University Press)
-  Haug H and Jauho A P 2008 Quantum Kinetics in Transport and Optics of Semiconductors (Springer)
-  Meir Y and Wingreen N S 1992 Phys. Rev. Lett. 68 2512
-  Dhar A and Shastry B S 2003 Phys. Rev. B 67 195405
-  Agarwal A and Sen D 2007 Phys. Rev. B 76 235316
-  Stefanucci G, Kurth S, Rubio A and Gross E K U 2008 Phys. Rev. B 77 075339
-  Tal-Ezer H and Kosloff R 1984 J. Chem. Phys. 81 3967
-  Alvermann A and Fehske H 2008 Phys. Rev. B 77 045125
-  Alvermann A and Fehske H 2009 Phys. Rev. Lett. 102 150601