Efficient auxiliary-mode approach for time-dependent nanoelectronics

The setup consists of a nanodevice coupled to two fermionic reservoirs. We adopt the usual partitioning of the total Hamiltonian into three terms: ${H}_{{\rm{T}}}(t)={H}_{\mathrm{dev}.}(t)+{H}_{\mathrm{res}.}(t)+{H}_{\mathrm{tun}.}$, where the device Hamiltonian is denoted by ${H}_{\mathrm{dev}.}$, the reservoirs (leads) are described by ${H}_{\mathrm{res}.}$, and ${H}_{\mathrm{tun}.}$ is the device-reservoir tunnel coupling term. More specifically, the device region is described by a general tight-binding (TB) Hamiltonian

$$\begin{eqnarray}&&{H}_{\mathrm{dev}.}(t)=\displaystyle \sum _{n=1}^{N}{\varepsilon }_{n}(t){c}_{n}^{\dagger }{c}_{n}-\displaystyle \sum _{n\ne m}^{N}{\gamma }_{{nm}}(t){c}_{n}^{\dagger }{c}_{m},\end{eqnarray} \tag{ 1 }$$

where ${c}_{n}^{\dagger }$ $({c}_{n})$ creates (annihilates) an electron at site n with on-site energy ${\varepsilon }_{n}$. Moreover, N is the total number of device sites and ${\gamma }_{{nm}}$ are the hopping elements between sites. The electron reservoirs consist of free particles and their Hamiltonian is given by

$$\begin{eqnarray}&&{H}_{\mathrm{res}.}(t)=\displaystyle \sum _{\alpha \in \{{\rm{L}},{\rm{R}}\}}\displaystyle \sum _{k}{\varepsilon }_{\alpha k}(t){b}_{\alpha k}^{\dagger }{b}_{\alpha k}.\end{eqnarray} \tag{ 2 }$$

Here, the subscript α denotes the left ($\alpha ={\rm{L}}$) or right ($\alpha ={\rm{R}}$) reservoir, k runs over the reservoirs states and ${b}_{\alpha k}^{(\dagger )}$ denote reservoir operators. The time-dependence of the reservoir energies ${\varepsilon }_{\alpha k}(t)$ is assumed to be given by ${\varepsilon }_{\alpha k}(t)={\varepsilon }_{\alpha k}^{0}+{{\rm{\Delta }}}_{\alpha }(t)$, i.e., all states in the reservoir are shifted by the same energy ${{\rm{\Delta }}}_{\alpha }(t)$. Lastly, the device-reservoir coupling is described by the Hamiltonian

$$\begin{eqnarray}&&{H}_{\mathrm{tun}.}=\displaystyle \sum _{n}\displaystyle \sum _{\alpha ,k}{T}_{{kn}}^{\alpha }{b}_{\alpha k}^{\dagger }{c}_{n}+{T}_{{kn}}^{\alpha * }{c}_{n}^{\dagger }{b}_{\alpha k},\end{eqnarray} \tag{ 3 }$$

where ${T}_{{kn}}^{\alpha }$ 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], ${T}_{{kn}}^{\alpha }(t)={T}_{{kn}}^{\alpha }{u}_{\alpha n}(t)$. For simplicity we consider ${T}_{{kn}}^{\alpha }$ as time-independent quantities in the following.

In general, the dynamics of the device coupled to electron reservoirs can be described by the single-particle RDO ${\boldsymbol{\sigma }}(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]

$$\begin{eqnarray}&&{\rm{i}}{\rm{\hslash }}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\boldsymbol{\sigma }}(t)=[{{\bf{H}}}_{\mathrm{dev}.}(t),{\boldsymbol{\sigma }}(t)]+{\rm{i}}\displaystyle \sum _{\alpha \in \{{\rm{L}},{\rm{R}}\}}({{\boldsymbol{\Pi }}}_{\alpha }(t)+{{\boldsymbol{\Pi }}}_{\alpha }^{\dagger }(t)),\end{eqnarray} \tag{ 4 }$$

where one defines the so-called current matrices ${{\boldsymbol{\Pi }}}_{\alpha }(t)$. Those can be expressed in terms of Green functions and self-energies

$$\begin{eqnarray}&&{{\boldsymbol{\Pi }}}_{\alpha }(t)={\int }_{{t}_{0}}^{t}{\rm{d}}{t}_{2}({{\bf{G}}}^{\gt }(t,{t}_{2}){{\boldsymbol{\Sigma }}}_{\alpha }^{\lt }({t}_{2},t)-{{\bf{G}}}^{\lt }(t,{t}_{2}){{\boldsymbol{\Sigma }}}_{\alpha }^{\gt }({t}_{2},t)),\end{eqnarray} \tag{ 5 }$$

where ${{\bf{G}}}^{\gtrless }$ and ${{\boldsymbol{\Sigma }}}_{\alpha }^{\gtrless }$ are the usual two-time greater/lesser Green's functions and self-energies [10], respectively. The self-energies are defined as [9]

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{\alpha }^{\gt }({t}_{2},t)=-{{\rm{ie}}}^{\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}{\int }_{{t}_{2}}^{t}{\rm{d}}{t}_{1}{{\rm{\Delta }}}_{\alpha }({t}_{1})}\int \displaystyle \frac{{\rm{d}}\epsilon }{2\pi {\rm{\hslash }}}(1-{f}_{\alpha }(\epsilon )){{\boldsymbol{\Gamma }}}_{\alpha }(\epsilon ){{\rm{e}}}^{-{\rm{i}}\epsilon ({t}_{2}-t)/{\rm{\hslash }}},\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{\alpha }^{\lt }({t}_{2},t)={{\rm{ie}}}^{\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}{\int }_{{t}_{2}}^{t}{\rm{d}}{t}_{1}{{\rm{\Delta }}}_{\alpha }({t}_{1})}\int \displaystyle \frac{{\rm{d}}\epsilon }{2\pi {\rm{\hslash }}}{f}_{\alpha }(\epsilon ){{\boldsymbol{\Gamma }}}_{\alpha }(\epsilon ){{\rm{e}}}^{-{\rm{i}}\epsilon ({t}_{2}-t)/{\rm{\hslash }}}.\end{eqnarray} \tag{ 6b }$$

Here, ${f}_{\alpha }(\varepsilon )\equiv f(\beta (\varepsilon -{\mu }_{\alpha }))$ is the Fermi distribution, which characterizes the equilibrium state of reservoir α with chemical potential ${\mu }_{\alpha }$ and inverse temperature $\beta ={({k}_{{\rm{B}}}T)}^{-1}$. Moreover, ${{\boldsymbol{\Gamma }}}_{\alpha }(\varepsilon )$ denotes the level-width function of reservoir α, which is defined in terms of the reservoir density of states ${\rho }_{\alpha }(\varepsilon )$

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\alpha {mn}}(\varepsilon )=2\pi {\rho }_{\alpha }(\varepsilon ){T}_{\alpha n}(\varepsilon ){T}_{\alpha m}^{* }(\varepsilon ).\end{eqnarray} \tag{ 7 }$$

The tunnel-coupling elements ${T}_{\alpha n}(\varepsilon )$ describe the coupling between device level n and the reservoir state at energy ε.

Once the current matrices, ${{\boldsymbol{\Pi }}}_{\alpha }(t)$ are known, the time-resolved current from reservoir α into the device is expressed in a very compact form [18]

$$\begin{eqnarray}&&{J}_{\alpha }(t)=\displaystyle \frac{2e}{{\rm{\hslash }}}\,\mathrm{Re}\,\mathrm{Tr}\{{{\boldsymbol{\Pi }}}_{\alpha }(t)\}.\end{eqnarray} \tag{ 8 }$$

Taking the trace of equation (4) and using the expression for J_α given above, yields the continuity equation for the electron current.

The auxiliary-mode expansion of ${{\boldsymbol{\Pi }}}_{\alpha }(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 ${{\boldsymbol{\Sigma }}}_{\alpha }^{\gtrless }$, 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.

2.3.1. 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

$$\begin{eqnarray}&&f(\varepsilon )=\displaystyle \frac{1}{1+\exp (\varepsilon )}\approx \displaystyle \frac{1}{2}-\displaystyle \sum _{p=1}^{{N}_{{\rm{F}}}}\left(\displaystyle \frac{{R}_{p}}{\varepsilon -{z}_{p}^{+}}+\displaystyle \frac{{R}_{p}}{\varepsilon -{z}_{p}^{-}}\right),\end{eqnarray} \tag{ 9 }$$

where R_p is the pth residue and ${z}_{p}^{\pm }$ is the pth pole in the upper $(+)$ and lower $(-)$ complex half-plane, whereas ${N}_{{\rm{F}}}$ is the total number of poles. Typical values for ${N}_{{\rm{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].

2.3.2. Decomposition of the level-width function

A central quantity for transport calculations is the level-width function ${{\boldsymbol{\Gamma }}}_{\alpha }(\varepsilon )$, introduced in equation (7). It can also be expressed in terms of the imaginary part of the retarded self-energy [10] via

$$\begin{eqnarray}&&{{\boldsymbol{\Gamma }}}_{\alpha }(\varepsilon )=-2\mathrm{Im}{{\boldsymbol{\Sigma }}}_{\alpha }^{{r}}(\varepsilon ).\end{eqnarray} \tag{ 10 }$$

The real part of ${{\boldsymbol{\Sigma }}}_{\alpha }^{{r}}$ is connected to ${{\boldsymbol{\Gamma }}}_{\alpha }$ via the Kramers–Kronig relation

$$\begin{eqnarray}&&\mathrm{Re}{{\boldsymbol{\Sigma }}}_{\alpha }^{{r}}(\varepsilon )=\displaystyle \frac{1}{\pi }{ \mathcal P }{\int }_{-\infty }^{\infty }{\rm{d}}\varepsilon ^{\prime} \displaystyle \frac{{{\boldsymbol{\Gamma }}}_{\alpha }(\varepsilon ^{\prime} )}{\varepsilon ^{\prime} -\varepsilon },\end{eqnarray} \tag{ 11 }$$

with ${ \mathcal P }$ 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 ${{\boldsymbol{\Sigma }}}_{\alpha }^{{r}}$ has to be computed numerically. Fortunately, many transport codes allow for the calculation of ${{\boldsymbol{\Sigma }}}^{{r}}$ [40].

Once the level-width function is known, it can be diagonalized [14, 39]

$$\begin{eqnarray}&&{{\boldsymbol{\Gamma }}}_{\alpha }(\varepsilon )=\displaystyle \sum _{c=1}^{{N}_{{\rm{c}}}}{\vec{\xi }}_{\alpha c}(\varepsilon ){v}_{\alpha c}(\varepsilon ){\vec{\xi }}_{\alpha c}^{\dagger }(\varepsilon ),\end{eqnarray} \tag{ 12 }$$

where ${\vec{\xi }}_{\alpha c}(\varepsilon )$ are transversal wave-vectors of dimension N, and ${v}_{\alpha c}(\varepsilon )$ 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 ${\vec{\xi }}_{\alpha c}(\varepsilon )$ 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, ${\vec{\xi }}_{\alpha c}$ 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 ${v}_{\alpha c}(\varepsilon )$ for each channel. Hence, the level-width function can be expressed by

$$\begin{eqnarray}&&{{\boldsymbol{\Gamma }}}_{\alpha }(\varepsilon )=\displaystyle \sum _{c=1}^{{N}_{c}}\,\displaystyle \sum _{{\ell }=1}^{{N}_{{\rm{L}}}}{\vec{\xi }}_{\alpha c}\displaystyle \frac{{{\rm{\Lambda }}}_{\alpha {\ell }c}{W}_{\alpha {\ell }c}^{2}}{{(\varepsilon -{\varepsilon }_{\alpha {\ell }c})}^{2}+{W}_{\alpha {\ell }c}^{2}}{\vec{\xi }}_{\alpha c}^{\dagger },\end{eqnarray} \tag{ 13 }$$

which is characterized by three (scalar) parameters ${{\rm{\Lambda }}}_{\alpha {\ell }c},{W}_{\alpha {\ell }c}$ and ${\varepsilon }_{\alpha {\ell }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 ${\vec{\xi }}_{\alpha c}$, one can also take advantage of the diagonalization. As in reference [18] the full level-width function is approximated by

$$\begin{eqnarray}&&{{\boldsymbol{\Gamma }}}_{\alpha }(\varepsilon )=\displaystyle \sum _{{\ell }=1}^{{N}_{{\rm{L}}}}\displaystyle \frac{{W}_{\alpha {\ell }}^{2}}{{(\varepsilon -{\varepsilon }_{\alpha {\ell }})}^{2}+{W}_{\alpha {\ell }}^{2}}{{\boldsymbol{\Lambda }}}_{\alpha {\ell }},\end{eqnarray} \tag{ 14 }$$

where ${{\boldsymbol{\Lambda }}}_{\alpha {\ell }}$ is now a matrix containing fitting coefficients. This matrix can be diagonalized and one obtains a set of eigenvectors and eigenvalues, for each α and ${\ell }$. The eigenvectors have an additional index ${\ell }$, but otherwise the procedure described in the next section also applies for this situation.

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

$$\begin{eqnarray}&&{{\boldsymbol{\Sigma }}}_{\alpha }^{\gtrless }({t}_{1},t)=\displaystyle \frac{1}{{\rm{\hslash }}}\displaystyle \sum _{c=1}^{{N}_{c}}\displaystyle \sum _{x=1}^{{N}_{{\rm{L}}}+{N}_{{\rm{F}}}}{\vec{\xi }}_{\alpha c}{{\rm{\Lambda }}}_{\alpha {xc}}^{\gtrless ,\pm }\exp \left\{\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}{\int }_{{t}_{1}}^{t}{\rm{d}}t^{\prime} ({\chi }_{\alpha {xc}}^{\pm }+{{\rm{\Delta }}}_{\alpha }(t^{\prime} ))\right\}{\vec{\xi }}_{\alpha c}^{\dagger }.\end{eqnarray} \tag{ 15 }$$

The new index x combines the Lorentz and Fermi pole indices, i.e., $x=\{{\ell },p\}$, where ${\ell }=1,\ldots ,{N}_{{\rm{L}}}$ and $p=1,\ldots ,{N}_{{\rm{F}}}$. Moreover, the coefficients entering (15) are

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\alpha {xc}}^{\gt ,\pm }=\left\{-\displaystyle \frac{{\rm{i}}}{2}{{\rm{\Lambda }}}_{\alpha {\ell }c}{W}_{\alpha {\ell }c}(1-{f}_{\alpha }({\varepsilon }_{\alpha {\ell }c}\pm {{\rm{i}}W}_{\alpha {\ell }c}))\,,\pm \displaystyle \frac{{R}_{p}}{\beta }{v}_{\alpha c}({\chi }_{\alpha {pc}}^{\pm })\right\},\end{eqnarray} \tag{ 16a }$$

$$\begin{eqnarray}&&{{\rm{\Lambda }}}_{\alpha {xc}}^{\lt ,\pm }=\left\{\displaystyle \frac{{\rm{i}}}{2}{{\rm{\Lambda }}}_{\alpha {\ell }c}{W}_{\alpha {\ell }c}{f}_{\alpha }({\varepsilon }_{\alpha {\ell }c}\pm {{\rm{i}}W}_{\alpha {\ell }c})\,,\pm \displaystyle \frac{{R}_{p}}{\beta }{v}_{\alpha c}({\chi }_{\alpha {pc}}^{\pm })\right\}\end{eqnarray} \tag{ 16b }$$

and

$$\begin{eqnarray}&&{\chi }_{\alpha {xc}}^{\pm }=\{{\varepsilon }_{\alpha {\ell }c}\pm {{\rm{i}}W}_{\alpha {\ell }c},{\mu }_{\alpha }+{z}_{p}^{\pm }/\beta \}.\end{eqnarray} \tag{ 17 }$$

Note that according to Jordan's lemma, if $t\gt {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\lt {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 ${{\boldsymbol{\Pi }}}_{\alpha }(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 wave-vectors

$$\begin{eqnarray}{\vec{{\rm{\Psi }}}}_{\alpha {xc}}(t) & = & {\displaystyle \int }_{{t}_{0}}^{t}{\rm{d}}{t}_{2}({{\bf{G}}}^{\gt }(t,{t}_{2}){{\rm{\Lambda }}}_{\alpha {xc}}^{\lt ,+}-{{\bf{G}}}^{\lt }(t,{t}_{2}){{\rm{\Lambda }}}_{\alpha {xc}}^{\gt ,+})\\ & & \times \exp \left\{\displaystyle \frac{{\rm{i}}}{{\rm{\hslash }}}{\displaystyle \int }_{{t}_{2}}^{t}{\rm{d}}t^{\prime} ({\chi }_{\alpha {xc}}^{+}+{{\rm{\Delta }}}_{\alpha }(t^{\prime} ))\right\}{\vec{\xi }}_{\alpha c},\end{eqnarray} \tag{ 18 }$$

which yields the decomposition of ${{\boldsymbol{\Pi }}}_{\alpha }(t)$ as follows

$$\begin{eqnarray}&&{{\boldsymbol{\Pi }}}_{\alpha }(t)=\displaystyle \frac{1}{{\rm{\hslash }}}\displaystyle \sum _{c=1}^{{N}_{c}}\displaystyle \sum _{x=1}^{{N}_{{\rm{L}}}+{N}_{{\rm{F}}}}{\vec{{\rm{\Psi }}}}_{\alpha {xc}}(t){\vec{\xi }}_{\alpha c}^{\dagger }.\end{eqnarray} \tag{ 19 }$$

By starting from their defining relation, equation (18), the following equation of motion for the auxiliary wave-vectors ${\vec{{\rm{\Psi }}}}_{\alpha {xc}}(t)$ is found

$$\begin{eqnarray}{\rm{i}}{\rm{\hslash }}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{\vec{{\rm{\Psi }}}}_{\alpha {xc}}(t) & = & [{{\bf{H}}}_{\mathrm{dev}.}(t)-({\chi }_{\alpha {xc}}^{+}+{{\rm{\Delta }}}_{\alpha }(t)){\bf{1}}]{\vec{{\rm{\Psi }}}}_{\alpha {xc}}(t)+{\rm{\hslash }}{{\rm{\Lambda }}}_{\alpha {xc}}^{\lt ,+}{\vec{\xi }}_{\alpha c}\\ & & +{\rm{\hslash }}({{\rm{\Lambda }}}_{\alpha {xc}}^{\gt ,+}-{{\rm{\Lambda }}}_{\alpha {xc}}^{\lt ,+}){\boldsymbol{\sigma }}(t){\vec{\xi }}_{\alpha c}+\displaystyle \frac{1}{{\rm{\hslash }}}\displaystyle \sum _{\alpha ^{\prime} x^{\prime} c^{\prime} }{{\rm{\Omega }}}_{\alpha {xc},\alpha ^{\prime} x^{\prime} c^{\prime} }(t){\vec{\xi }}_{\alpha ^{\prime} c^{\prime} },\end{eqnarray} \tag{ 20 }$$

where terms containing time derivatives of the greater and lesser Green's functions are contained in the two-mode quantity ${{\rm{\Omega }}}_{\alpha {xc},\alpha ^{\prime} x^{\prime} c^{\prime} }(t)$. Its definition is given in the appendix. The initial condition of (20), i.e. ${\vec{{\rm{\Psi }}}}_{\alpha {xc}}({t}_{0})=0$, follows from equation (18). In a similar manner, the time evolution of ${{\rm{\Omega }}}_{\alpha {xc},\alpha ^{\prime} x^{\prime} c^{\prime} }(t)$ is obtained as

$$\begin{eqnarray}{\rm{i}}{\rm{\hslash }}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{{\rm{\Omega }}}_{\alpha {xc},\alpha ^{\prime} x^{\prime} c^{\prime} }(t) & = & ({\chi }_{\alpha ^{\prime} x^{\prime} c^{\prime} }^{-}(t)+{{\rm{\Delta }}}_{\alpha ^{\prime} }(t)-{\chi }_{\alpha {xc}}^{+}(t)-{{\rm{\Delta }}}_{\alpha }(t)){{\rm{\Omega }}}_{\alpha {xc},\alpha ^{\prime} x^{\prime} c^{\prime} }(t)\\ & & +{\rm{i}}{\rm{\hslash }}({{\rm{\Lambda }}}_{\alpha ^{\prime} x^{\prime} c^{\prime} }^{\gt ,-}-{{\rm{\Lambda }}}_{\alpha ^{\prime} x^{\prime} c^{\prime} }^{\lt ,-}){\vec{\xi }}_{\alpha ^{\prime} c^{\prime} }^{\dagger }{\vec{{\rm{\Psi }}}}_{\alpha {xc}}(t)+{\rm{i}}{\rm{\hslash }}({{\rm{\Lambda }}}_{\alpha {xc}}^{\gt ,+}-{{\rm{\Lambda }}}_{\alpha {xc}}^{\lt ,+}){\vec{{\rm{\Psi }}}}_{\alpha ^{\prime} x^{\prime} c^{\prime} }^{\dagger }(t){\vec{\xi }}_{\alpha c},\end{eqnarray} \tag{ 21 }$$

with ${{\rm{\Omega }}}_{\alpha {xc},\alpha ^{\prime} x^{\prime} c^{\prime} }({t}_{0})=0$. It is crucial to note that in the present approach the quantities ${{\rm{\Omega }}}_{\alpha {xc},\alpha ^{\prime} x^{\prime} c^{\prime} }$ 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. Moreover, ${{\rm{\Omega }}}_{\alpha {xc},\alpha ^{\prime} x^{\prime} c^{\prime} }(t)=0$ if $x={N}_{{\rm{L}}}+1,\ldots ,{N}_{{\rm{L}}}+{N}_{{\rm{F}}}$ and $x^{\prime} ={N}_{{\rm{L}}}+1,\ldots ,{N}_{{\rm{L}}}+{N}_{{\rm{F}}}$, that is if x and $x^{\prime}$ represent auxiliary modes resulting from the Fermi-function expansion. This is easily understood by noting that the last two additive terms in (21) vanish since ${{\rm{\Lambda }}}_{\alpha ^{\prime} x^{\prime} c^{\prime} }^{\gt ,\pm }={{\rm{\Lambda }}}_{\alpha ^{\prime} x^{\prime} c^{\prime} }^{\lt ,\pm }$ for those $x,x^{\prime}$ named above. Then (21) is directly integrable yielding ${{\rm{\Omega }}}_{\alpha {xc},\alpha ^{\prime} x^{\prime} c^{\prime} }(t)\propto {{\rm{\Omega }}}_{\alpha {xc},\alpha ^{\prime} x^{\prime} c^{\prime} }({t}_{0})=0$, 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 (4), (19)–(21) constitute a closed system of ODEs, which can be solved by standard methods. The device Hamiltonian ${{\bf{H}}}_{\mathrm{dev}.}$ and the energy shift in the leads ${{\rm{\Delta }}}_{\alpha }$ can have an arbitrary time-dependence. By propagating the equations one obtains the RDO, and via equation (8), the currents between device and leads at each time step. In total, there are ${N}_{\mathrm{leads}}\cdot {N}_{{\rm{c}}}\cdot ({N}_{{\rm{L}}}+{N}_{{\rm{F}}})$ auxiliary wave-vectors of length N and ${N}_{\mathrm{leads}}^{2}\cdot {N}_{{\rm{c}}}^{2}\cdot {N}_{{\rm{L}}}\cdot ({N}_{{\rm{L}}}+2{N}_{{\rm{F}}})$ scalars to be propagated. For comparison, in the original approach [18], there are ${N}_{\mathrm{leads}}\cdot ({N}_{{\rm{L}}}+{N}_{{\rm{F}}})$ auxiliary current matrices (N × N) and ${N}_{\mathrm{leads}}^{2}\cdot {N}_{{\rm{L}}}\cdot ({N}_{{\rm{L}}}+2{N}_{{\rm{F}}})$ ${\boldsymbol{\Omega }}$-matrices (N × N). The computational complexity is therefore drastically reduced in the present wave-vector method.

Our first numerical investigation deals with the multi-Lorentzian decomposition of $v(\varepsilon )$. In figure 1 we consider the effect of truncating the approximation (13) at a certain upper bound ${N}_{{\rm{L}}}$. The $3\times {N}_{{\rm{L}}}$ parameters ${{\rm{\Lambda }}}_{\alpha {\ell }},{\varepsilon }_{\alpha {\ell }},{W}_{\alpha {\ell }}$ are determined by a least squares fit to expression (23), where $v(\varepsilon )$ is evaluated at discrete energies. The resulting fit is shown in figure 1(b) together with the exact expression. One can see that in particular at the band edges $(\varepsilon =\pm 2\gamma )$ the approximated level-width function displays oscillations and attains negative values. However, with increasing number of Lorentzians these artifacts become less pronounced. Already for a rather small number of Lorentzians the agreement is very good (see, e.g., the red curve displaying the case ${N}_{{\rm{L}}}=16$).

To assess the effect of the decomposition of the level-width function on the overall calculation, we show in figure 1(c) a current–voltage curve. In effect, we calculate the steady-state current for different ${N}_{{\rm{L}}}$ and varying voltage. The steady-state current is obtained by letting the system evolve to the equilibrium state. Specifically, the simulation time amounts to ${t}_{\mathrm{sim}}=1500\,{\rm{\hslash }}/\gamma$. Note that in this example the voltage drop is applied by symmetrically shifting the energy bands of the leads, i.e., ${{\rm{\Delta }}}_{{\rm{L}}}={eV}/2$ and ${{\rm{\Delta }}}_{{\rm{R}}}=-{eV}/2$, at equal chemical potentials (${\mu }_{{\rm{L}}}={\mu }_{{\rm{R}}}=0$). 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., ${eV}={\mu }_{{\rm{L}}}-{\mu }_{{\rm{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}_{{\rm{L}}}=16$.

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}_{{\rm{F}}}=60$, to have a strong effect of memory load, rather than for reasons of convergence (${N}_{{\rm{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 ${ \mathcal O }({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.

Zoom In Zoom Out Reset image size

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 ${ \mathcal O }({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].

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 time-dependent 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, ${\mu }_{{\rm{L}}}={\mu }_{{\rm{R}}}=0$. Starting from equilibrium, a Lorentzian voltage pulse, ${{\rm{\Delta }}}_{{\rm{L}}}{(t)=2{\rm{\hslash }}\tau /((t-{t}_{0})}^{2}+{\tau }^{2})$, with pulse length $\tau =7.5{\rm{\hslash }}/\gamma$ is applied to the left lead. The time t₀ 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 $\mu \mathrm{eV}$. 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, $({I}_{{\rm{L}}}-{I}_{{\rm{R}}})/2$, 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., ${q}_{\alpha }\equiv {\int }_{-\infty }^{t}{I}_{\alpha }(t^{\prime} ){\rm{d}}t^{\prime}$, is equal to the injected or extracted charge. We have checked that ${q}_{{\rm{L}}}\approx -{q}_{{\rm{R}}}\approx 1$ for each of the above cases. For small system sizes ($N\lesssim 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.

Zoom In Zoom Out Reset image size

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}_{i}=1/2$. 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 ${v}_{\mathrm{Lev}}\approx 2\gamma a/{\rm{\hslash }}$. This is in accordance with the semi-classical picture of the electrons, which gives for the group velocity [47]

$$\begin{eqnarray}&&{v}_{\mathrm{Lev}}=\displaystyle \frac{1}{{\rm{\hslash }}}\displaystyle \frac{\partial }{\partial k}\varepsilon (k){| }_{\varepsilon =\mu }=\displaystyle \frac{2\gamma a}{{\rm{\hslash }}}\,.\end{eqnarray} \tag{ 24 }$$

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].

Zoom In Zoom Out Reset image size