Efficient auxiliary-mode approach for time-dependent nanoelectronics

A new scheme for numerically solving the equations arising in the time-dependent non-equilibrium Green's function formalism is developed. It is based on an auxiliary-mode expansion of the self-energies which convert the complicated set of integro-differential equations into a set of ordinary differential equations. In the new scheme all auxiliary matrices are replaced by vectors or scalars. This drastically reduces the computational effort and memory requirements of the method, rendering it applicable to topical problems in electron quantum optics and molecular electronics. As an illustrative example we consider the dynamics of a Leviton wave-packet in a 1D wire.


Introduction
The new field of electron quantum optics [1][2][3] attracts a lot of interest and increases the demand for the theoretical description of time-resolved phenomena in nanoscale devices. Recent experimental developments allow time-dependent investigations and manipulation of electronic excitations above the Fermi sea, which had been predicted by Levitov and coworkers [4,5]. Understanding the influence of external or environmental driving on the transport and dynamics of electrons is also important for other systems, such as molecular junctions [6] or nanoelectromechanical systems [7].
However, theoretical descriptions of such situations remain very challenging. Several methods are available, but efficient approaches for general-purpose computations are still under investigation. A popular approach in this regard is the (time-dependent) nonequilibrium Greenʼs function (TDNEGF) formalism [8][9][10], which provides an exact treatment of the electron transport under time-dependent conditions. Since the resulting equations are hard to solve numerically, there are many attempts for dealing with them. In particular, several 'brute-force' approaches have been reported [11][12][13][14] as well as a recursive technique [15]. There exist also other theories, which are mathematically equivalent to TDNEGF, e.g., the partition-free approaches developed in [16,17], or the wave-function-based scheme of reference [14].
In [18], the authors use an auxiliary-mode expansion to recast the TDNEGF theory in a tractable form. This scheme is based on two main ingredients: a multi-Lorentzian parametrization [19] of the level-width function as well as the Matsubara or an equivalent decomposition of the Fermi function [20][21][22]. In essence, these expansions in terms of simple poles result in a decomposition of the self-energies as a sum of weighted exponentials. For each term in the sum, an auxiliary current-matrix can be defined and the equations of motion of those are easily obtained. Finally, a closed system of ordinary differential equations (ODEs) for the singleparticle reduced density operator (RDO) and the auxiliary quantities is derived. In this formalism, memoryeffects are completely accounted for and there is no restriction on the magnitude of the device-reservoir coupling. Nevertheless, the computational time scales linearly with the number of time-steps. This method, along with its close variations (see, e.g., [23] for the general case and [24] for the wide-band limit), has already been successfully applied to study several realistic scenarios [25][26][27][28][29] and to model systems [30,31].
The bottleneck of the auxiliary-mode method results from dealing with many auxiliary N×N matrices, where N is the number of sites in the system. The resulting memory and computational requirements prevent the application of the method to larger systems. In this article we present a new version of the auxiliary-mode propagation scheme, where all auxiliary matrices are replaced by N-dimensional vectors or scalars. This is achieved by diagonalizing the level-width function and using the Lorentzian decomposition for the available transport channels. As a consequence, the computational effort and memory requirements of the method are drastically reduced. It opens the possibility to study a variety of systems, from molecular junctions to electron quantum optics.
The article is organized as follows. In section 2 we derive the new vector-based method. Specifically, in section 2.1 we describe the transport setup and the Hamiltonian. Sections 2.2 and 2.3 provide details of the auxiliary-mode approach, on which the present method is based. Lastly, in section 2.4 we discuss the new formulation of the theory in terms of wave-vectors and show the final form of the equations of motion. Section 3 is dedicated to numerical results. First, we investigate the multi-Lorentzian decomposition in section 3.1 and in section 3.2 we show a performance analysis of the new method. In section 3.3 the scheme is applied to the propagation of a Leviton wave-packet in a 1D system. Finally, section 4 gives a summary and conclusions.

Hamiltonian
The setup consists of a nanodevice coupled to two fermionic reservoirs. We adopt the usual partitioning of the total Hamiltonian into three terms: tun. , where the device Hamiltonian is denoted by H dev. , the reservoirs (leads) are described by H res. , and H tun. is the device-reservoir tunnel coupling term. More specifically, the device region is described by a general tight-binding (TB) Hamiltonian where † c n ( ) c n creates (annihilates) an electron at site n with on-site energy e n . Moreover, N is the total number of device sites and g nm are the hopping elements between sites. The electron reservoirs consist of free particles and their Hamiltonian is given by Here, the subscript α denotes the left (a = L) or right (a = R) reservoir, k runs over the reservoirs states and where a T kn are tunneling elements, which connect device state n to reservoir state k. Note that in principle the tunneling elements can be time-dependent as well. This is achieved by imposing a factorized form [9,18] For simplicity we consider a T kn as time-independent quantities in the following.

Reduced density operator
In general, the dynamics of the device coupled to electron reservoirs can be described by the single-particle RDO s ( ) t . Here and in the following, bold-face variables denote matrices in the device state-space which is spanned by N basis functions. The time-evolution of the RDO is given by an equation of the following form [32] å s s where one defines the so-called current matrices P a ( ) t . Those can be expressed in terms of Green functions and self-energies where  G and S a  are the usual two-time greater/lesser Greenʼs functions and self-energies [10], respectively.
The self-energies are defined as [9] ò ò p Here, The tunnel-coupling elements e a ( ) T n describe the coupling between device level n and the reservoir state at energy ε.
Once the current matrices, P a ( ) t are known, the time-resolved current from reservoir α into the device is expressed in a very compact form [18] Taking the trace of equation (4) and using the expression for J α given above, yields the continuity equation for the electron current.

Auxiliary-mode expansion
The auxiliary-mode expansion of P a ( ) t introduced in reference [18] is based on a twofold decomposition: first, the Fermi function is expanded in terms of simple poles in the complex plane and secondly the level-width function is expressed as a sum of Lorentzians. The main goal of the auxiliary-mode expansion is the conversion of integro-differential equations, like equation (4), into a system of ODEs. This is achieved by approximating the memory kernel in equation (5), which is given by the self-energies S a  , by a sum of exponential functions. Such schemes have successfully been used in the context of non-Markovian quantum master equations [19,33] or within the hierarchical equations of motion theory [34,35]. Moreover, reference [36] proposes also other classes of functions to be used for this parametrization in the case when the self-energies exhibit super-ohmic behavior.

Decomposition of the Fermi function
There exist several possibilities to expand the Fermi function into a sum of simple poles [20,21,37]. Best known is the Matsubara decomposition [37], which has, however, very poor convergence properties and is not practical for low temperatures. A very efficient expansion is the Padé decomposition [20], which leads to å e e e e = + » -- where R p is the pth residue and  z p is the pth pole in the upper + ( )and lower -( )complex half-plane, whereas N F is the total number of poles. Typical values for N F used in the present article range from 20 to 60 poles. Note that the residues and poles can be efficiently calculated by solving an eigenvalue problem [38].

Decomposition of the level-width function
A central quantity for transport calculations is the level-width function e G a ( ), introduced in equation (7). It can also be expressed in terms of the imaginary part of the retarded self-energy [10] via with  denoting the Cauchy principal value. The calculation of the retarded self-energy for a given system-lead configuration is in general non-trivial. For some cases analytical expressions are available [39], but in principle S a r has to be computed numerically. Fortunately, many transport codes allow for the calculation of S r [40].
Once the level-width function is known, it can be diagonalized [14,39] å and e a ( ) v c is directly proportional to the group velocity in channel c. The factor of proportionality is a/ÿ with a the typical lattice constant in the leads. Moreover, N c denotes the total number of transport channels. For transport geometries the number of transport channels is typically much lower than the dimension N of the problem, as it is approximately given by the number of links between the lead and the system. The eigenvectors x e a ( ) c  in (12) are in general energy dependent, normalized and non-orthogonal [14]. However, in certain cases, such as a square-lattice (or simple cubic) structure in the leads without magnetic field, x ac  become independent of energy [39]. Those situations are of particular importance and we focus on energy independent wave-vectors in the next sections, but note that the scheme also works in the general case. It is crucial to note that for a given geometry the diagonalization (12) is done only once at the beginning of the simulation.
While in reference [18] the full level-width function is approximated by a sum of Lorentzians, in the present case we only have to expand e a ( ) v c for each channel. Hence, the level-width function can be expressed by W , c c and e aℓc . These parameters can be obtained by a partitioned nonlinear regression [41]. The combination of the Lorentzian decomposition and the wave-vector expansion is the key for the improved scheme, which leads to a drastic simplification of the method.
In the more general case of energy-dependent eigenvectors x ac  , one can also take advantage of the diagonalization. As in reference [18] the full level-width function is approximated by where L aℓ is now a matrix containing fitting coefficients. This matrix can be diagonalized and one obtains a set of eigenvectors and eigenvalues, for each α and ℓ. The eigenvectors have an additional index ℓ, but otherwise the procedure described in the next section also applies for this situation.

Auxiliary-mode wave-vectors
As a next step, we use the decomposition of the level-width and Fermi functions, equations (13) and (9), and insert them into the definitions of the greater/lesser self-energies, i.e., (6a) and (6b). The energy integration can be performed employing Jordanʼs lemma, which brings the greater/lesser self-energies into the following form The new index x combines the Lorentz and Fermi pole indices, i.e., = ℓ Note that according to Jordanʼs lemma, if > t t 1 then those coefficients above bearing the index + must be used. This stands for poles (and their corresponding residues) from the upper half plane. Conversely, when < t t 1 one uses those parameters bearing the index −.
As one can see, the self-energies are given by a sum of weighted exponentials. Now we use the expansion (15) in the definition of the current matrices P a ( ) t , which is given by equation (5). Note that the Greenʼs functions under the integral are standing to the left of the self-energies. This allows us to define auxiliary-mode wavevectors  is obtained as . It is crucial to note that in the present approach the quantities W a a¢ ¢ ¢ xc x c , are scalars, rather than matrices as in the original scheme [18]. This renders the method very efficient in terms of computational time and memory consumption.
F , that is if x and ¢ x represent auxiliary modes resulting from the Fermi-function expansion. This is easily understood by noting that the last two additive terms in (21) , in accordance with our initial condition. This observation is also consistent with the fact, that the retarded self-energy does not depend on the Fermi function. Equations

Results and discussion
In order to demonstrate the proposed method, we discuss in the following numerical results for a 1D-TB chain. In this case the Hamiltonian introduced in section 2.1 takes the following form Note that we consider the same hopping mechanism for the device as well as for the reservoirs, with positionindependent hopping constant γ. This means that partitioning between system and reservoirs regions occurs only by correspondingly indexing the sites (see figure 1(a)). For the 1D case there is one transport channel per lead ( = N 1 c ), which means that the level-width function has only one non-zero element, given by e ( ) v [39] e g e g e g For the following examples we use an adaptive Runge-Kutta method [42] to solve the ODEs. Then, from the resulting RDO we extract the occupation numbers as s = n i i i , and from the auxiliary wave-vectors we obtain the currents via equation (8). ). This situation is termed as a purely electric voltage drop. It is to be distinguished from the case when the voltage is applied as a difference in chemical potentials, i.e., m m = -eV L R , termed as purely chemical voltage drop [9,14]. Since there is no on-site potential in the device, it is sufficient to consider only a few sites for the current calculation. The device simulated in figure 1(c) consists of two sites. The resulting stationary current shows a deviation with respect to the exact result in the range of large voltages. Those deviations arise due to the aforementioned artefacts found for the fitted level-width function. Accordingly, they disappear on this scale when using = N 16 L .

Numerical efficiency
Next, we explore the reduction of computation time achieved by the present method. The run-times are compared to the approach derived in [18]. Specifically, in figure 2 the scaling of the computational time with system size N is shown, for N up to 256 sites. Note that we choose a large number of Padé terms, i.e., = N 60 F , to have a strong effect of memory load, rather than for reasons of convergence (N F as low as 20 suffices for simulations with the given parameters). The data clearly show the superior performance of the newly developed scheme. Although both methods eventually scale as ( )  N 3 , the reduction in computational effort by the new scheme is evident. This results from the speed-up and lower memory requirement of matrix-vector multiplications as compared to matrix-matrix multiplications in the original method.
The limiting factor on the route to an even more extensive reduction of computational effort is given by equation (4), which is present in both methods. The latter contains the commutator of the RDO with the Hamiltonian, it thus involves matrix-matrix multiplications. For a larger system size N this operation becomes computationally expensive, which is seen in figure 2. For small N the computation time scales approximately linearly with N (blue curve). In contrast, for larger N it approaches the expected ( )  N 3 dependence, which results from dense matrix-matrix multiplications. On the other hand, having the density matrix available at each step of the computation is one main advantage of the present method compared to other schemes, like the pure wave-vector method of [14].

Propagation of leviton wave-packets
Finally, as an application of our method, we discuss the propagation of Leviton wave-packets through a 1D-TB chain. Levitons are collective many-body states, which involve the excitation of electrons above the Fermi sea, with no corresponding creation of particle-hole pairs. Such minimal-excitation states, theoretically predicted by Levitov and coworkers [4,5] and recently measured experimentally (see, e.g., [3]), are observed when a timedependent potential is applied to the system initially in equilibrium. In the particular case of Lorentzian-shaped voltage-pulses, it turns out that the associated current-noise vanishes, which indicates the absence of additional hole excitations. Such noiseless single electron sources have been proposed for the realization of electron quantum optics [3,43,44].
In the simulation we consider 1D-TB chains, with up to 128 sites, coupled to left and right leads at zero chemical potentials, m m = = 0 L R . Starting from equilibrium, a Lorentzian voltage pulse, is applied to the left lead. The time t 0 is taken as a reference time ( = t 0 0 ) in the following. This pulse generates a coherent electronic wave-packet carrying exactly one elementary charge [4].
A note on the physical time scales of such a scenario: Lorentzian pulses in the range of tens of picoseconds were reported in [3], which facilitated the generation of Leviton excitations in a quasi-1D system. Specifically, the studied device was a quantum point contact defined by lateral gates applied to a two-dimensional electron gas (2DEG) [3]. Such pulses correspond in our model to hopping elements γ on the order of 50 meV. In addition, [45] proposes using a 2DEG in the quantum Hall regime to study Levitons, whereas in [46] the authors theoretically discuss the generation of Levitons in graphene.
In figure 3 we show the time-resolved net current, - , through our device for various system sizes ( = N 16,32,64,128). Corresponding to the Leviton entering as well as exiting the device, one expects peaks in the left and right currents, respectively. Moreover, the area under those peaks, i.e., is equal to the injected or extracted charge. We have checked that » -» q q 1 L R for each of the above cases. For small system sizes (  N 16), the two events showing the wave-packet entering and leaving the system cannot be discerned from each other, due to the finite temporal width of the pulse. In this case one obtains approximately one peak in the net current (see the lowest dashed curve in figure 3). By increasing the system size the temporal spacing between the two peaks becomes larger, permitting to distinguish the two events. It is this case that we want to focus on in the following. Figures 4(a) and (c) show the left and right time-resolved currents, respectively, for N=128. In the time between entering and leaving, the Leviton crosses the system and can therefore be manipulated. The actual propagation through the device is visualized by the time-dependent occupancy of the sites relative to their equilibrium value, = n 1 2 i . This is displayed in figure 4(b), which shows the trace of the Leviton through the device. Eventually, the negative peak in the right current, in figure 4(c), indicates that the Leviton is leaving the device and entering the right lead. One can see that the center of the wave-packet follows a straight trajectory through the device. From the time difference between the current peaks we estimate the group velocity of the Leviton to be . This is in accordance with the semi-classical picture of the electrons, which gives for the group velocity [47] This illustrative example shows that the present method is able to easily treat a system which is large enough (N=128) and to investigate the Leviton propagation within a reasonable computation time. Such a scenario is not feasible with the method of [18].

Conclusions
We proposed a new scheme for the propagation of the reduced density matrix of a nanoelectronic system. Specifically, by diagonalizing the level-width function of the leads, the complicated integro-differential equations of the TDNEGF method can be transformed into ODEs for wave-vectors and scalars. Employing the latter quantities instead of the large-size auxiliary matrices encountered previously results in a significant reduction of memory load. This further sets the stage for a drastic decrease of computation time, which allows it to treat larger systems. We demonstrated the method by showing numerical results for one-dimensional systems. However, the new scheme is perfectly suited to treat also other geometries, such as two-dimensional square or hexagonal lattices. In order to achieve good performance, the number of transport channels has to be smaller than the number of sites in the device. Or, in other words, the level-width function has to have only a few non-zero eigenvalues. This is indeed the case in many transport situations, since the number of sites coupled to leads is typically much smaller than the size of the system. Asymptotically, the scaling of the method is linear in the number of time-steps and cubic in the system-size. The latter can be further reduced if the device Hamiltonian only contains short-range couplings. The cubic scaling is the main bottleneck for treating much larger systems. On the other hand, the RDO and the currents are available at each time-step without extra cost. This allows it to investigate scenarios where the device Hamiltonian depends on the RDO. For example, electron-electron interactions on the mean-field level can easily be accounted for. It is also an advantage in the context of mixed quantum mechanical/molecular mechanical simulations, where the forces acting on the atoms depend on the electronic state of the device [48,49].
As an example, we treated the propagation of a Leviton wave-packet through a 1D-TB chain composed of 128 sites. Such a scenario, where a quasi-particle excitation well defined in time-domain is investigated, is becoming increasingly interesting as recent experiments start to focus on time-resolved quantum nanoelectronics [2,3]. Our method allows it to investigate the influence of interactions and de-phasing on the Leviton propagation. For example, the latter may be realized by using appropriate stochastic fluctuations of the site energies [50]. Overall, our scheme considerably extends the range of applications of the auxiliary-mode approach and offers a versatile method to investigate time-resolved phenomena in nanoelectronic devices.