ORDER O(1) ALGORITHM FOR FIRST-PRINCIPLES CALCULATION OF TRANSIENT CURRENT THROUGH OPEN QUANTUM SYSTEMS

Abstract

A fast algorithm is used to study the transient behavior due to the step-like pulse. This algorithm consists of two parts: The algorithm I reduces the computational complexity to T0N3 for large systems as long as T

Description

FIELD OF THE INVENTION

This invention relates to a transient current calculation as part of the study of the response time and peak transient current in the design of nano-chips in nano electronics.

BACKGROUND OF THE INVENTION

At the heart of growing demands for nanotechnology is the need of ultrafast transistors whose response time is one of the key performance indicators. The response of a general quantum open system can be probed by sending a step-like pulse across the system and monitoring its transient current over time, thus making transient dynamics a very important consideration. Many experimental data show that most molecular device characteristics are closely related to the material and chemical details of the device structure. Therefore, first principles analysis, which involves quantitative and predictive analysis of device characteristics, especially their dynamic properties, without relying on any phenomenological parameter, becomes a central interest in the design of nanoelectronics.

The theoretical study of transient current dates back twenty years to when the exact solution in the wideband limit (“WBL”) was obtained by Wingreen et al. See, N. S. Wingreen et al., Phys. Rev. B 48, 8487(R) (1993), which is incorporated herein by reference in its entirety. Since then transient current has been studied extensively using various methods, including the scattering wave function, non-equilibrium Green's function (“NEGF”) approach, and density matrix method. See, J. Wang, J. Comput. Electron. 12, 343-355 (2013); S. Kurth et al., Phys. Rev. B 72, 035308 (2005)(hereinafter “Kurth”); G. Stefanucci et al, Phys. Rev. B 77, 075339 (2008)(hereinafter “Stefanucci”); Y. Zhu et al., Phys. Rev. B 71, 075317 (2005); J. Maciejko et al., Phys. Rev. B 74, 085324 (2006)(hereinafter “Maciejko”); R. Tuovinen et al., Phys. Rev. B 89, 085131 (2014)(herein after “Tuovinen”); R. Seoane Souto et al., Phys. Rev. B 92, 125435 (2015) and X. Zheng et al., Phys. Rev. B 75, 195127 (2007), all of which are incorporated herein by reference. The major obstacle to the theoretical investigation of transient current is its computational complexity. Many attempts have been made trying to speed up the calculation See Kurth and Tuovinen, as well as L. Zhang et al., Phys. Rev. B 86, 155438 (2012); B. Gaury et al., Phys. Rep. 534, 1-37 (2014); M. Ridley et al., J. Phys.: Conf. Ser. 696, 012017 (2016)(hereinafter “Ridley”); H A. Croy et al., Phys. Rev. B 80, 245311 (2009) and F J. Weston et al., Phys. Rev. B 93, 134506 (2016), all of which are incorporated herein by reference in their entirety.

Despite these efforts, the best algorithm for calculating the transient current from first principles going beyond WBL limit scales like TN³ using a complex absorbing potential (“CAP”), where T and N are the number of time steps and the size of the system, respectively. See, L. Zhang et al, Phys. Rev. B 87, 205401 (2013) (hereinafter “Zhang”), which is incorporated herein by reference in its entirety. It should be noted that if WBL is used, the scaling is reduced. See, the Ridley article. However, to capture the feature of the band structure of the lead and the interaction between the lead and the scattering region, WBL is not a good approximation in the first principles calculation. As a result, most of the first principles investigations of transient dynamics were limited to small and simple one-dimensional systems.

There are a number of problems, such as with magnetic tunneling junctions (MTJ) and ferroelectric tunneling junctions, where the system is two dimensional or even three dimensional in nature. See, Z. Y. Ning et al., Phys. Rev. Lett. 100, 056803 (2008) and J. D. Burton et al., Phys. Rev. Lett. 106, 157203 (2011), which are incorporated herein by reference in their entirety. For these systems, a large number of k points N_k have to be sampled in the first Brillouin to capture accurately the band structure of the system. For MTJ structures like Fe—MgO—Fe, at least N_k=10⁴ k points must be used to give a converged transmission coefficient. See, D. Waldron et al., Phys. Rev. Lett. 97, 226802 (2006), which is incorporated herein by reference in its entirety. This makes the time consuming transient calculation N_k times longer, which is an almost impossible task even with high performance supercomputers. Clearly it is urgent to develop better algorithms to reduce the computational complexity.

SUMMARY OF THE INVENTION

According to the present invention a novel algorithm based on NEGF-CAP formalism is used to calculate transient current as a function of time step T. The computational time of this algorithm is independent of T as long as T<N² where N is the system size. Hence the algorithm is order O(1) as long as T<N². Four important factors are essential to achieve this result: (1) the availability of an exact solution of transient current based on a non-equilibrium Green's function (NEGF) that goes beyond the wideband limit (WBL); (2) the use of complex absorbing potential (CAP) so that the transient current can be expressed in terms of poles of Green's function; (3) within the NEGF-CAP formalism the transient current can be calculated separately in the space and time domains making the O(1) algorithm possible so that at this point the computational complexity reduces to 50N³+TN² (algorithm I); and (4) the exploitation of a Vandermonde matrix enables the use of the fast multipole method (FMM) and fast Fourier transform (FFT) to further reduce the scaling to 50N³+2N² log₂ N for T<N² and large N, which is therefore completely independent of T (algorithm II). See V. Rokhlin, J. Comput. Phys., 60, 187-207 (1985)(hereinafter “Rokhlin”) and J. Song et al., IEEE trans. Antennas Propagat. 45, 1488-1493 (1997)(hereinafter “Song”), which are both incorporated herein by reference in their entirety. It is possible with this new algorithm to calculate, and thus predict, the actual response of a nano-device theoretically without experimental data. Thus, the algorithm is useful in the design of a nano-device before actual production.

Everything is made of atoms. In theory, every physical phenomena in a real device can be accurately calculated and predicted once its composition is known, e.g., how many atoms, which types of atoms, etc. From the Hamiltonian (a quantity that stores the physical information of the system based on the particular electrons of particular elements present in the devices), it is possible to use the algorithm of the present invention to calculate and predict the actual response of a nano-device. To verify the computational complexity, benchmark calculations were carried out on graphene nanoribbons using the tight-binding model. A speed up factor of 1000 T was gained for a system size of N=2400. A calculation was also done for the same system with N=10,200 and T=10⁸ confirming the O(1) scaling. This fast algorithm makes the computational complexity of transient current calculation comparable to that of a static calculation. Thus, the huge speed gain enables the performance of first principles transient calculations and nano-device design on a modest computer workstation.

For a general open quantum system with multiple leads under a step-like bias pulse, the Hamiltonian (the operator corresponding to the total energy of the system in most of the cases) is given by

$H\left(t\right)=\sum _{k\alpha }{\epsilon }_{k\alpha }{\hat{c}}_{k\alpha }^{†}{\hat{c}}_{k\alpha }+\sum _n\left({\epsilon }_n+U_n\left(t\right)\right){\hat{d}}_n^{†}{\hat{d}}_n+\sum _{k\alpha n}h_{k\alpha n}{\hat{c}}_{k\alpha }^{†}{\hat{d}}_n+c.c.,$

where c^† (c) denotes the electron creation (annihilation) operator in the lead region. The first term in this equation corresponds to the Hamiltonian of leads with ε_kα the energy of lead α which contains external bias voltage v_α(t)=V_αθ(±t). The second and third terms represent the Hamiltonian in the central scattering region and its coupling to leads, respectively. Also included is the time-dependent Coulomb interaction U_n(t) in the scattering region. This Hamiltonian can be obtained using the first-principle method or it can be assumed to be a tight-binding form. Taking q=h=1, the time-dependent terminal current I_α(t) of lead a is defined as in the Zhang article.

I _α(t)=ReTr[┌_α(2H(t)−i∂_t)G^<(t,t)┌_α],  (1)

where ┌_α is an auxiliary projection matrix which is used for measuring the transient current passing through the lead α. Here G^< and H(t) are the lesser Green's function and the Hamiltonian of the central scattering region, respectively. For a step-like pulse an exact solution for G^< has been obtained by Maciejko et al. as described in the Maciejko article, which goes beyond the WBL and has been applied to first principles calculation of transient current for atomic junctions. Next, the case of an upward pulse v_α(t)=V_αΘ(t) can be considered. In order to have an exact solution as described in the Maciejko article, it is assumed that U_n(t)=U_n,eq+(U_n,neq−U_n,eq)Θ(t) where U_n,eq is the equilibrium Coulomb potential, while U_n,neq is the non-equilibrium potential at a long time limit. Note that the fast algorithm relies on this approximation where the self-consistent Coulomb potentials are needed only at two points: t=0⁻ and t=∞. As a result of this approximation, at t=0, the system is in equilibrium with Hamiltonian H_eq while at t>0 the system is in a non-equilibrium state with a time-independent Hamiltonian H_neq. In terms of a spectral function A_α(ε,t), the lesser Green's function G^<is given by the Maciejko article as

$\begin{array}{cc}{G}^{<}\left(t,t\right)=i\sum _{\alpha }\int ·d\varepsilon 2\pi f\left(\varepsilon \right)A_{\alpha }\left(\varepsilon ,t\right){\Gamma }_{\alpha }\left(\varepsilon \right)A_{\alpha }^{†}\left(\varepsilon ,t\right).& \left(2\right)\end{array}$

The upward step-like bias pulse the A_α(ε,t) is found in the Maciejko article to be

$\begin{array}{cc}\begin{array}{c}{A}_{\alpha }\left(\epsilon ,t\right)={\stackrel{\_}{G}}^r\left(\epsilon +{\Delta }_{\alpha }\right)-\int \frac{d\omega }{2\pi i}\frac{e^{-i\left(\omega -\epsilon \right)t}{\stackrel{\_}{G}}^r\left(\omega +{\Delta }_{\alpha }\right)}{\left(\omega -\epsilon +{\Delta }_{\alpha }-i0^{+}\right)}\times \\ \left[\frac{{\Delta }_{\alpha }}{\left(\omega -\epsilon -i0^{+}\right)}+\Delta {\tilde{G}}^r\left(\epsilon \right)\right]\\ \equiv A_{1\alpha }\left(\epsilon +{\Delta }_a\right)+\int d\omega e^{-i\left(\omega -\epsilon \right)t}A_{2\alpha }\left(\omega ,\epsilon \right),\end{array}& \left(3\right)\end{array}$

where G^r and {tilde over (G)}^r are the non-equilibrium and equilibrium retarded Green functions, respectively, Δ_α is the amplitude of the external bias −V_α, and Δ=U_neq−U_eq is a matrix where the subscript “neq” and “eq” refer to non-equilibrium and equilibrium potentials, respectively.

Despite the simplification from the conventional double time G^<(t,t′) to single time G^<(t,t) used in Eq. (1Error! Reference source not found.), the computational cost to obtain G^< remains very demanding for the following reasons:

(1) Consider A_α(ε,t) with a matrix size of N, matrix multiplications G^r(ω+Δ_α) and {tilde over (G)}(ε)^r in the integrand of Eq. (3) require computational complexity of O(N³) for each time step. As a result, the total computational cost over a period of time is at least O(TN³) where T is the number of time steps.

(2) Double integrations in energy space are required for G^<. The presence of numerous quasi-resonant states whose energies are close to the real energy axis makes the energy integration in A_α extremely difficult to converge.

This problem can be overcome using the CAP method as described in the article J. Driscoll and K. Varga, Phys. Rev. B 78, 245118 (2008), which is incorporated herein by reference in its entirety. The essence of the CAP method is to replace each semi-infinite lead by a finite region of CAP while keeping the transmission coefficient of the system unchanged. In addition, it is demonstrated in the Zhang article that the first principles result of transient current for molecular junctions obtained from the exact numerical method (non-WBL) and the CAP method are exactly the same. Using the CAP method, the poles of Green's function can be obtained easily and the spectral function can be calculated analytically using the residue theorem. Expanding the Fermi function using the Pade spectrum decomposition (PSD), as explained more fully in the detailed description of the invention below, further allows for the calculation of the transient current separately in space and time domain making the O(T⁰N³) algorithm possible. See, J. Hu et al., J. Chem. Phys. 133, 101106 (2010), which is incorporated herein by reference in its entirety.

Achieving the algorithm of the present invention for the transient current calculation can now be demonstrated, i.e., I_α(t_j) for j=1,2, . . . ,T. First Eq. (3) is substituted into Eq.(2), and G^<(t,t) can be written as

$\begin{array}{cc}{G}^{<}\left(t,t\right)=\left(i/\pi \right)\left[B_1+\int d\omega d{\omega }^{\prime }e^{-i\left(\omega -{\omega }^{\prime }\right)t}B_2\left(\omega ,{\omega }^{\prime }\right)+\sum _{\alpha }\int d\varepsilon d{\omega }^{\prime }e^{i\left({\omega }^{\prime }-\epsilon \right)t}f\left(\epsilon \right)A_{1\alpha }W_{\alpha }A_{2\alpha }^{†}+c.c.\right],& \left(4\right)\end{array}$

where,

B₁=∫dεf(ε)Σ_αA_1αW_αA_1α^†, B₂(ω,ω^t)=∫dεf(ε)Σ_αA_2α(ω,ε)W_αA_2α^†(ω^t,ε), and W_α is the CAP matrix. In terms of the poles of Green's function and the Fermi distribution function, which are explained more fully in the detailed description of the invention below, we have:

$\begin{array}{cc}{G}^{<}\left(t,t\right)=\left(i/\pi \right)\left[B_1+\sum _{\mathrm{nm}}e^{-i\left({\epsilon }_n-{\epsilon }_m^{*}\right)t}{\stackrel{\_}{B}}_2\left({\epsilon }_n,{\epsilon }_m^{*}\right)++\sum _{\alpha }\sum _{\mathrm{lm}}e^{-i\left({\epsilon }_l-{\epsilon }_m^{*}+{\Delta }_{\alpha }\right)t}\stackrel{\_}{f}\left({\stackrel{\_}{\epsilon }}_l\right){\stackrel{\_}{B}}_{3\alpha }\left({\stackrel{\_}{\epsilon }}_l,{\epsilon }_m^{*}\right)+c.c.+\sum _{\alpha }\sum _{\mathrm{nm}}e^{-i\left({\epsilon }_n-{\epsilon }_m^{*}\right)t}f\left({\epsilon }_n-{\Delta }_{\alpha }\right){\stackrel{\_}{B}}_{4\alpha }\left({\epsilon }_n,{\epsilon }_m^{*}\right)+c.c.\right],& \left(5\right)\end{array}$

where ε_n and ε_m (n=1,2, . . . N) is the complex energy spectrum of

$H_{\mathrm{neq}}-i\sum _{\alpha }W_{\alpha }$

in the lower half plane, while {tilde over (ε)}_l are the poles of f(E) using PSD with l=1, . . . N_f: and N_f is the total number of those poles for the adopted Pade approximant.

Within the CAP framework, G^< in Eq. (1) is the lesser Green's function of the central scattering region excluding the CAP regions. Substituting the second term of Eq. (4) into the first term in Eq. (1), its contribution to current (denoted as I₁) is

$\begin{array}{cc}\begin{array}{c}{I}_1\left(t\right)=\mathrm{Re}\sum _{\mathrm{nm}}e^{-i\left({\epsilon }_n-{\epsilon }_m^{*}\right)t}\mathrm{Tr}\left[{\stackrel{\_}{\Gamma }}_{\alpha }2H_{\mathrm{neq}}{\stackrel{\_}{B}}_2\left({\epsilon }_n,{\epsilon }_m^{*}\right){\stackrel{\_}{\Gamma }}_{\alpha }\right]\\ \equiv 2\mathrm{Re}\sum _{\mathrm{nm}}e^{-i\left({\epsilon }_n-{\epsilon }_m^{*}\right)t}M_{\mathrm{nm}},\end{array}& \left(6\right)\end{array}$

where the matrix M does not depend on time. Thus, the space and time domains have been separated.

Denoting a Vandermonde matrix V_jk=exp(iε_kt_j) with k=1,2, . . . ,N, t_j=jdt, j=1,2, . . . ,T where dt is the time interval and I₁(t_j)=[V^t(M+M^†)V*]_jj. Using this approach, I_L is finally obtained. This calculation is more fully explained in the detailed description below.

$\begin{array}{cc}{I}_L\left(t_j\right)=I_{0L}+{\left[V^tM_1V^{*}\right]}_{\mathrm{jj}}+\left(\sum _{\alpha =L,R}{\left[V^tM_{2\alpha }{\stackrel{\_}{V}}_{\alpha }^{*}\right]}_{\mathrm{jj}}+c.c.\right),& \left(7\right)\end{array}$

where {circumflex over (V)}_αjk=exp(i({circumflex over (ε)}_k+Δ_α)t_j) is a T×N_f matrix, M₁ is a N×N matrix while M₂ is a N×N_f matrix. Since ε_k is the complex energy in the lower half plane, V_jk goes to zero at large j. Hence I_0L is the long time limit of the transient current which can be calculated using the Landauer Buttiker formula. The time dependent part of the transient current can be separated into a real space calculation (calculation of M₁ and M_2α) and then a matrix multiplication involving time. At room temperatures the Fermi function can be accurately approximated by 15 or 20 Pade approximants. Hence the calculation of (V^tM₁V*|_ff+(Σ_α|V^tM_2α{tilde over (V)}_α*|_ff+α,α) can be combined to give TN² computational complexity.

The computational complexity of this algorithm (denoted as algorithm I) can now be examined By counting the number of matrix multiplication, the computational complexity of the real space calculations can be estimated to be 50N³. Therefore the total computational complexity is 50N³+TN². At this stage, the algorithm is not O(1) yet. In the detailed description below it is shown that matrix multiplication V^tM, where M is M₁ or M_2α, can be done using the FMM and FFT (denoted as algorithm II). This reduces the computational complexity of V^tM from TN² to T log₂N. Hence for T<N², the computational complexity is 50N³+N² log₂N. For T>N², the scaling is 50N³+T log₂N. However, for large T, the physics comes into play. Since ε_j is the complex energy of the resonant state, V_jT=exp(−iε_jdtT) decays quickly to zero before T=N². For a graphene nanoribbon with N=10⁴ (see details below), the maximum value of V_jT=exp(−iε_jdtT) is 10⁻³ when T=N and dt=1 fs. Consequently all of the matrix elements are zero for T=10N. Hence for large systems, the chance to go beyond T=N² is small. In this sense, algorithm II is an order O(1) algorithm.

BRIEF DESCRIPTION OF THE DRAWINGS

The foregoing and other objects and advantages of the present invention will become more apparent when considered in connection with the following detailed description and appended drawings in which like designations denote like elements in the various views, and wherein:

FIG. 1A shows the configuration of a gated graphene nanoribbon and a voltage applied to it over time and FIG. 1B shows the transient current of the zigzag graphene nanoribbon for a system of 600 atoms for the present invention and the Zhang method;

FIG. 2 shows the scaling of computation time versus the number of computations N for the prior art Zhang method and the Algorithm I of the method according to the present invention;

FIG. 3 shows the scaling of computational time against time t for N equal to 10,200, time step dt=1 fs;

FIG. 4 shows a graph of the typical behavior of transient current, wherein the inset is the long time behavior between 17 ps to 18 ps;

FIG. 5 is a graph of the transmission coefficient of a zigzag graphene nanoribbon for a system of 10,000 atoms calculated by 30 layers of CAP and by exact numerical solution;

FIG. 6 shows the distribution of (1/ω_T)^j on the complex plane;

FIG. 7 shows the distribution of d_j on the complex plane;

FIG. 8 illustrates the transient current calculated by Algorithm I, Algorithm IIa and Algorithm IIb; and

FIG. 9 is a flow chart showing the process for the calculation of the algorithm according to the present invention.

DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTS

First principles transient current calculation is essential to the study of the response time and to capture the peak transient current for preventing the melt down of nano-chips in nano-electronics. For a period of time T, its calculation is known to be extremely time consuming with the best scaling being TN³ where N is the dimension of the device. The dynamical response of the system is usually probed by sending a step-like pulse and monitoring its transient behavior.

The present invention is directed to a new algorithm for calculating the transient current of a system, which requires much less computational time than the prior art. FIG. 9 shows a flow chart of the method for calculating the algorithm of the present invention. It can be carried out on a computer workstation as part of the design of a nano-device.

Generally speaking, the algorithm provides an O(1) computational method for obtaining the transient response of current I(t) over the entire user-defined time period under suitable conditions. The input parameters for this algorithm include Hamiltonians of the open quantum system before and after the transient process denoted as H_eq and H_neq. These Hamiltonians, whether based on atomistic first-principle density functional theory (DFT) or tight-binding methods, can be adopted at will for the system of interest. Nonetheless, complex absorbing potential (CAP) is required to represent the self-energy of the leads of the device. Subsequently, corresponding eigenstate and energy can be generated and fed into the algorithm. Additionally, Fermi function is expanded using Pade spectrum decomposition (PSD).

It is only possible to construct the particular exact equation of transient current which can be calculated separated in space and time domains by using the mentioned form of input ingredients. Once such an equation is constructed, whether for the case for step-up, step-down bias or other similar classes of transient response, so that I(t) is a function of I_o, V_m(t) and M_n for some m and n where V_m(t) are Vandemonde matrices and I_o and M_n are time independent matrices. Once the required Vandemonde matrices and M_n are constructed, fast multi pole method (FMM) and fast Fourier transform (FFT) are used to calculate the matrix multiplication between them. Eventually, I(t) can be obtained by summing up the contribution due to I_o, V_m(t) and M_n.

In step 901 of FIG. 9 the process begins with input preparation. This includes step 901A, which involves forming Preliminary Input quantities H_eq, H_neq (From first-principle (DFT), TB, etc) and Self energy of leads (from CAP). In step 901B eigenstates and eigen-energy are generated. This input information is then directed to step 902. In step 902 construction is performed on the exact equation of transient current based on NEGF formalism that is organized into the general form of I(t)=I_t(I_o, V_m(t), M_n) such that it is a function composed of time dependent Vandermonde matrices V_m and other space dependent matrices M_n. Also, there is an input from step 903, i.e., Fermi functions (from the PSD method) to step 902.

The results of step 902 are then reviewed in step 904 to see if they contain time dependent components. If the answer is NO, the process moves to step 905. In step 905 the process computes and constructs space dependent matrices by an optimized matrix multiplication process. From step 905 the process goes to step 906 where the FMM and FFT methods are used for computing the multiplication of V_m and M_n.

If the result of step 904 is YES, the process moves to step 907 where Vandermonde matrices are prepared before moving to step 906. The results in both step 905 and 906 are summed in Step 908. The summed results are directed to output step 909, where the complete transient response I(t) over the user-defined time period is generated.

To demonstrate the power of the algorithm according to the present invention, the transient current in a graphene nanoribbon is calculated. Graphene is a well-known intrinsic 2D material with many exotic properties. See, A. H. Castro Neto et al., Rev. Mod. Phys. 81, 109 (2009) and Y. Zhang et al., Nature, 438, 201-204 (2005), which are incorporated herein by reference in their entirety. Studies of its transient behavior in response to a step-like pulse have been reported in the literature. See the Stefanucci article and E. Perfetto et al., Phys. Rev. B 82, 035446 (2010) and Y. O. Klymenko et al., Eur. Phys. J. B 69, 383-388 (2009), which are both incorporated herein by reference in their entirety.

The algorithm of the present invention was tested on a gated graphene nanoribbon at room temperature using the tight-binding (TB) Hamiltonian given by

$\begin{array}{cc}\hat{H}=-h_0\sum _{<i,j>}{\hat{c}}_i^{†}{\hat{c}}_j-q\sum _i\left[V_1\theta \left(t\right)+V_{g1i}+V_{g2i}\right]{\hat{c}}_i^{†}{\hat{c}}_j,& \left(8\right)\end{array}$

where ĉ_i^†(ĉ_i) is the creation (annihilation) operator at site i and h₀=2.7 eV being the nearest hopping constant. Here V(x)=V_L+(V_R−V_L)×/L is the potential landscape due to the external bias with V_R=−V_L=0.54V and V_g1 and V_g2 being the gate voltages in regions S₁ and S₂, respectively.

First it is confirmed that the transient current calculated using the new method is the same as that of the Zhang reference. Using 30 layers of CAP, the transmission coefficient versus the energy was calculated and it showed good agreement with the exact solution. FIG. 5 provides a comparison for a graphene nanoribbon with N=10,000. In FIG. 5 the solid line (T_S) is the exact numerical result using self-energy of the lead and the dashed line (T_CAP) is from CAP. This also ensures the correct steady state current. For the transient current, excellent agreement is also obtained between the algorithm of the present invention and that of the Zhang reference (See FIG. 2). It should be noted that in the presence of gates, the “on-off” time of graphene is shortened in comparison to un-gated graphene which has a long oscillating current as reported in the Stefannucci article. Also the transient current is calculated for a graphene nanoribbon with size of N=10⁴ and T=20,000 with each time step ifs. See FIG. 5.

Now testing of the scaling of the algorithm of the present invention can be achieved by calculating the transient current for nanoribbons with different system sizes ranging from 600 to 10,200 atoms. The first test is with algorithm I. The computational time for the transient current for 3 time steps compared against system sizes N is shown in FIG. 2. In FIG. 5 the data is fitted using 50N³+TN² with very good agreement showing TN² scaling for the time-dependent part. For comparison, FIG. 5 also shows a plot of the computation time using the method in the Zhang reference. It was found that the number of energy points N_E depends on the spectrum of resonant states of the system. For graphene nanoribbons with 600 atoms, N_E=6000 was used to converge the integral over a Fermi function. FIG. 2 shows that a speed up factor of 1000T is achieved at N=2400. The scaling is shown in FIG. 3, from which it can be seen that for T<N the computational time is almost independent of the number of time steps.

Now algorithm II, which reduces the scaling TN² further, can be examined Notice that the scaling TN² comes from matrix multiplication involving Vandermonde matrix V^tM₁. The fast algorithm is available to speed up the calculation involving a structured matrix, such as the Vandermonde matrix. As discussed in detail below the FMM as in the Rokhlin and Song articles and FFT can be used to carry out the same matrix multiplication using only c₃N² log₂N operations provided that T<N². Here the coefficient c₃ is a large constant that depends only on the tolerance of the calculation τ and the setup of FMM. The theoretical estimate of this coefficient is about 40 log₂(1/τ) where τ is the tolerance in the FMM calculation, which used τ=10⁻⁴. See N. Yarvin et al., Anal. 36, 629 (1999), which is incorporated herein by reference in its entirety. When implementing FMM, this coefficient is in general larger than the theoretical one.

To test algorithm II, the transient current is calculated for N=10⁴ and T=10⁸ as explained in detail below, using FMM and FFT. Denote t₁ the CPU time needed for the spacial calculation (50N³), t₂ is the time needed for the temporal part (matrix multiplication in Eq. (7)) using, e.g., a Xeon X5650 workstation with 12 cores and a frequency of 2.67 GHz. The result t₁=3500 s is obtained using 12 cores and t₁=33800 s using a single core so the efficiency of multithreading is about 80%. For an FMM calculation, multithreading could be very inefficient so a single core has been used to perform the calculation. It was found t₂=3400 s for T=10⁸ using a single core, as shown below in the detailed computational analysis and numerical calculation. It was found that for T=10⁸ the time spent in the time dependent part is about one tenth of the time of the independent calculation. This confirms that the method of the present invention uses an algorithm of order O(1) as long as T<N². It should be pointed out that algorithm II is directed to the calculation of the transient current l(t) with time steps T=N² at one shot with scaling N² log₂N. This scaling remains if l(t) with the number of time step less than N² is desired.

FIG. 2 shows the scaling of N against computation time at T=3 (each time step is 1 fs). The fitted curve, which is in the form of 50N³+TN², is in good agreement with the calculated results The Y-axis is on the right. In order compare the present invention to the performance reported in the Zhang reference, 6000 energy points were used for integration using the Y-axis on the left.

FIG. 3 illustrates the scaling of t against CPU time for N=10,200 (t=100 fs corresponds to T=100) using algorithm I. The left hand side uses an exponential scale in t and the right hand side uses a linear scale in t. FIG. 3 shows that at extremely large data points, e.g., ranges over 10 thousands points, the computational time is proportional to T.

Since the algorithm of the present invention is based on the NEGF-CAP formalism, it can be extended to the NEGF-DFT-CAP formalism which performs the first principles calculation. In fact, the NEGF-DFT-CAP method has already been successfully implemented in the first principles transient current calculation as shown in the Zhang reference, which gives exactly the same result as that of NEGF-DFT. With the fast algorithm at hand, many applications can be envisaged. For instance, the transient spin current (related to spin transfer torque) using the NEGF-DFT-CAP formalism has been carried out for planar structures where k-sampling in the first Brillouin zone is needed. It is straightforward to include k-sampling in the method of the present invention. It is also possible to extend this method to the case where electron-phonon interaction in the Born approximation as well as other dephasing mechanisms are present. Finally, first principles transient photo-induced current on two dimensional layered materials can be calculated using the method of the present invention.

Some of the details of the calculations presented above are given here.

Pade Approximant

Brute force integration over the Fermi function along the real energy axis to obtain G^<(t,t) may need thousands of energy points to converge, which is very inefficient. To obtain an accurate result while reducing the cost, fast converging Pade spectrum decomposition (PSD) is used for the Fermi function f in Eq. (4) above, so that the residue theorem can be applied. Using [n−1/n] PSD scheme with the Pade approximant accurate up to O((ε/kT)⁴ⁿ⁻¹), Fermi function f can be expressed as

$\begin{array}{cc}f\left(\varepsilon \right)=\frac{1}{2}-\sum _{j=1}^n\frac{2{\eta }_j\mathrm{\beta \varepsilon }}{{\left(\mathrm{\beta \varepsilon }\right)}^2+{\xi }_j^2},& \left(9\right)\end{array}$

where and ζ_j and n_j are two set of constants that can be calculated easily. Using the PSD scheme analytic form of G^< in Eq. (4) can be obtained using the residue theorem. See J. Hu et al., J. Chem. Phys. 133, 101106 (2010), which is incorporated herein by reference in its entirety.

Calculation of the Spectral Function

The terms and Ĝ^r(ε) and G^r(ε), the equilibrium and non-equilibrium retarded

Green's functions, respectively, can be expressed in terms of their eigen-functions by solving the following eigen-equations for H_eq and H_neq, i.e.,

(H_eq−iW)ψ_n⁰=ε_n⁰ψ_n⁰,

(H_eq+iW^†)φ_n⁰=ε_n⁰φ_n⁰,  (10)

where

$W=\sum _{\alpha }W_{\alpha }$

and similar equations can be defined for H_neq. See the Zhang reference. Using the eigen-functions of H_eq−iW and H_neq−iW, we have

$\begin{array}{cc}{\tilde{G}}^r\left(\epsilon \right)={\left[\epsilon -H_{\mathrm{eq}}+\mathrm{iW}\right]}^{-1}=\sum _n\frac{{\psi }_n^0〉〈{\phi }_n^0}{\left(\epsilon -{\epsilon }_n^0+i0^{+}\right)},& \left(11\right)\\ {\stackrel{\_}{G}}^r\left(\epsilon \right)={\left[\epsilon -H_{\mathrm{neq}}+\mathrm{iW}\right]}^{-1}=\sum _n\frac{{\psi }_n〉〈{\phi }_n}{\left(\epsilon -{\epsilon }_n+i0^{+}\right)}.& \left(12\right)\end{array}$

Performing an integral over ω using the residue theorem, the analytic solution of A_α is obtained

$\begin{array}{cc}{A}_{\alpha }\left(\epsilon ,t\right)=\sum _n\frac{{\psi }_n〉〈{\phi }_n}{\epsilon +{\Delta }_{\alpha }-{\epsilon }_n+i0^{+}}+\sum _n\frac{e^{{i\left(\epsilon +{\Delta }_{\alpha }-{\epsilon }_n\right)}^t}{\psi }_n〉〈{\phi }_n}{\epsilon -{\epsilon }_n+i0^{+}}\times \hspace{1em}\left[\frac{{\Delta }_{\alpha }}{\epsilon +{\Delta }_{\alpha }-{\epsilon }_n+i0^{+}}-\Delta \sum _l\frac{{\psi }_l^0〉〈{\phi }_l^0}{\left(\epsilon -{\epsilon }_l^0+i0^{+}\right)}\right],& \left(13\right)\end{array}$

where Δ=H_neq−H_eq.

Calculation of the Lesser Green's Function

In Eq. (5) using residue theorem the involved terms are defined as

${\stackrel{\_}{B}}_2=-4{\pi }^2\left[B_2\left(\omega ,{\omega }^{\prime }\right)\left(\omega -{\epsilon }_n\right)\left({\omega }^{\prime }-{\epsilon }_m^{*}\right)\right]{}_{\omega ={\epsilon }_n,{\omega }^{\prime }={\epsilon }_m^{*}}{\stackrel{\_}{B}}_{3\alpha }=-2\pi {\mathrm{iA}}_{1\alpha }\left({\tilde{\epsilon }}_l\right)W_{\alpha }\left[A_{2\alpha }^{†}\left({\omega }^{\prime },{\tilde{\epsilon }}_l\right)\left({\omega }^{\prime }-{\epsilon }_m^{*}\right)\right]{}_{{\omega }^{\prime }={\epsilon }_m^{*}}{\stackrel{\_}{B}}_{4\alpha }=-4{\pi }^2[A_{1\alpha }\left(\epsilon \right)W_{\alpha }A_{2\alpha }^{†}\left({\omega }^{\prime },\epsilon \right)\left({\omega }^{\prime }-{\epsilon }_m^{*}\right)\times \hspace{1em}\left(\epsilon -{\epsilon }_n+{\Delta }_{\alpha }\right)]{}_{\epsilon ={\epsilon }_n-{\Delta }_{\alpha },{\omega }^{\prime }={\epsilon }_m^{*}}\stackrel{\_}{f}=2\pi i\left(f\left(\epsilon \right)\left(\epsilon -{\tilde{\epsilon }}_l\right)\right){}_{\epsilon ={\stackrel{\_}{\epsilon }}_l}$

Calculation of the Transient Current

Starting from Eq. (1) and in analogue to Eq. (6), the expressions of the current in Eq. (7) can be obtained as follows:

$I_{0L}\left(t_j\right)=2\mathrm{Re}\mathrm{Tr}\left[\frac{i}{\pi }{\stackrel{\_}{\Gamma }}_LH_{\mathrm{neq}}B_1{\stackrel{\_}{\Gamma }}_L\right]$ $M_1=\mathrm{Re}\mathrm{Tr}\left[\frac{i}{\pi }{\stackrel{\_}{\Gamma }}_L\left(2H_{\mathrm{neq}}-\left({\epsilon }_n-{\epsilon }_m^{*}\right)\right)\left({\stackrel{\_}{B}}_2+\sum _{\alpha }f\left({\epsilon }_n-{\Delta }_{\alpha }\right){\stackrel{\_}{B}}_{4\alpha }\right){\stackrel{\_}{\Gamma }}_L\right]$ $M_{2\alpha }=\mathrm{Re}\mathrm{Tr}\left[\frac{i}{\pi }{\stackrel{\_}{\Gamma }}_L\left(2H_{\mathrm{neq}}-\left({\stackrel{\_}{\epsilon }}_l-{\epsilon }_m^{*}+{\Delta }_o\right)\right)\left(f\left({\epsilon }_m^{*}\right){\stackrel{\_}{B}}_{3\alpha }\right){\stackrel{\_}{\Gamma }}_L\right].$

The expression of transient current I_R(t) is similar to Eq. (7).

Transient Current for a Araphene Nanoribbon of System Size N=10⁴

A calculation was performed on transient current through a zigzag graphene nanoribbon of 10,000 atoms with T=20,000 time steps (each time step is 1 fs). The width of the system is two unit cells (16 atoms) while the length of the system is 625 unit cells. Two gate voltages of 2.2V were applied so that the system is in the tunneling regime. The bias voltage is v_L=−v_R=0.5 V. From FIG. 4, it can be seen that a typical behavior of transient current with the current shooting up initially and then decreasing oscillations to the long time limit (the numerical results using FMM show that 100 ps is needed to reach the dc limit, which is shown as a dashed line). The oscillatory behavior is due to resonant states in the system. This test was at a temperature of 300° K. The inset is the long time behavior between 17 ps to 18 ps.

FIG. 5 shows the transmission coefficient of the zigzag graphene nanoribbon for a system of 10,000 atoms. The solid line (T_S) is the exact numerical result using self-energy of the lead and the dashed line (T_CAP) is from CAP. Thus, FIG. 5 shows a comparison between the exact numerical result of the transmission coefficient using the self-energy of the lead and that of CAP. There is very close agreement.

Fast Multipole Method

The fast multipole method has been widely used and has been ranked top 10 best algorithms in 20th Century. See, V. Rokhlin, J. Comput. Phys. 60, 187-207 (1985); J. Song, C. C. Lu, and W. C. Chew, IEEE trans. Antennas Propagat. 45, 1488-1493 (1997); and B. A. Cipra, SIAM News, 33(4), 2 (2000), which are incorporated herein by reference in their entirety. It is extremely efficient for large N. The following quantity is then calculated:

$\begin{array}{cc}I\left(t\right)=\sum _{n,m}\exp \left(-i{\varepsilon }_nt\right)M_{\mathrm{nm}}\exp \left(i{\varepsilon }_m^{*}t\right),& \left(14\right)\end{array}$

where the matrix M can be expressed in terms of vectors as M=(c₀,c₁, . . . ,c_N−1) and V_nj=exp(−iε_nt_j) is a Vandermonde matrix with t_j=jdt and j=1,2, . . . T. Eq. (14) is of the form V^tMV* where t stands for transpose. In the following, an outline of how to calculate V^tc where c is a vector of N components is given.

Setting a_j=exp(−iε_jdt) and denoting T the number of time steps. Then b=V^tc is equivalent to

$b_n=\sum _{j=0}^{N-1}{c_j\left(a_j\right)}^n.$

A direct computation shows that the entries of b=V^tc are the first T coefficients of the Taylor expansion of

$\begin{array}{cc}S\left(x\right)=\sum _{j=0}^{N-1}\frac{c_j}{1-a_jx}=\sum _n^{\infty }\sum _{j=0}^{N-1}{c_j\left(a_jx\right)}^n=\sum _nb_nx^n,& \left(15\right)\end{array}$

where b_n=Σ_j=0^N−1 c_j(a_j)ⁿ. Denoting S(x)=Σ_m=0^T−1 b_nx^m and setting x=ω_T^l with ω_T=exp(i2π/T) it can be used to calculate ⁻S(ω_T^l) which is the Fourier transform of b_n,

$\begin{array}{c}S\left({\omega }_T^l\right)=\sum _{j=0}^{N-1}\sum _{n=0}^{T-1}c_ja_j^n{\omega }_T^{\mathrm{nl}}\\ =\sum _{j=0}^{N-1}c_j\frac{1-{\left(a_j{\omega }_T^l\right)}^T}{1-a_j{\omega }_T^l}\\ ={\omega }_T^{-l}\sum _{j=0}^{N-1}\frac{c_j\left(1-a_j^T\right)}{{\left(1/{\omega }_T\right)}^l-a_j},\end{array}$

where ω_T^T=1 is used. Note that the fast multipole method (FMM) aims to calculate

$v_l=\sum _jc_j/\left(x_l-a_j\right)$

with O(N) operations instead of N² operations. Hence S(ω_T^l) can be obtained using FMM, from which b_n can be calculated using FFT.

Now the computational complexity for T≦N can be estimated. For FMM the value κ₁max(T,N) operations are needed where κ₁ is about 40 log₂(1/τ) with τ the tolerance. See N. Yarvin and V. Rokhlin, SIAM J. Numer. Anal. 36, 629 (1999), which is incorporated herein by reference in its entirety. For FFT the computational complexity is at most κ₂N log₂N, where κ₂ is a coefficient for FFT calculation. To compute V^tM where M has N vectors, V^tc is calculated N times. Hence the total computational complexity is κ₁N²+κ₂N² log₂N. This algorithm is denoted as algorithm IIa while the algorithm for T<N² discussed below is denoted as algorithm IIb.

For very large T up to T=N² (if N=10⁴ and T=10⁸), it can be shown that the computational complexity is κ₁N²+2κ₂N² log₂N. In fact, it is easy to see that I(t_j) defined in Eq. (6) is the first T coefficients of the Taylor expansion of

$\begin{array}{cc}\begin{array}{c}S\left(x\right)=\sum _{n,m=0}^{N-1}\frac{M_{\mathrm{nm}}}{1-a_na_m^{*}x}\\ =\sum _j^{\infty }\sum _{n,m=0}^{N-1}{M_{\mathrm{nm}}\left(a_na_m^{*}\right)}^ix^j\left(17\right)\\ =\sum _jI\left(t_j\right);\end{array}& \left(16\right)\end{array}$

where a_n=exp(−ε_ndt). Now two new vectors u and d can be defined which have N² components with u^t=(c₀^t,c₁^t, . . . ,c_N−1^t) (recall the definition M=(c₀,c₁, . . . ,c_N−1)) and d^t=(a₀*a^t,a₁*a^t, . . . ,a_N−1*a^t), where once again t stands for transpose. With the new vectors defined, S(x) in Eq. (16) is expressed as

$\begin{array}{cc}S\left(x\right)=\sum _{j=0}^{N^2-1}\frac{u_j}{1-d_jx},& \left(18\right)\end{array}$

which is exactly the same form as Eq. (15). The only difference is that c and a in Eq. (15) have N components and S has to be calculated N times while u and d in Eq. (18) have N² components and S can be calculated according to Eq. (18) just once. Therefore the computational complexity is κ₁N²+κ₂N² log₂N². If T=nN with n=1,2, . . . N, it is not difficult to show that the computational complexity is κ₁TN/n+κ₂T(N/n) log₂(nN)=κ₁N²+κ₂N² log₂(nN).

To summarize, the computational complexity of Eq. (14) is κ₁N²+2κ₂N² log₂N for T<N². It is easy to show that for T>N² the scaling is κ₁N²+2κ₂T log₂N. However, for large T, the physics comes into play. Since a_j=exp(−ε_jdt) with ε_j the energy of resonant state, a_j^T quickly decays to zero before T=N² and hence no need to go up for T>N².

Algorithm II was tested numerically for a system with N=10⁴ and T=10⁸. The configuration of the system is the same as that which appears in FIG. 1, except the width W of the system is now 17 times wider with a total of 10200˜10⁴ atoms. Each time step is 0.012 fs. The computed transient current using Algorithm I and II are shown in FIG. 8. The purpose of this calculation is to test the computational complexity only. All that needs to be done is to compute,

$\begin{array}{cc}{\hspace{0.17em}}^{-}S\left({\omega }_T^l\right)={\omega }_T^{-l}\sum _{j=0}^{N^2-1}\frac{u_j\left(1-d_j^T\right)}{{\left(1/{\omega }_T\right)}^l-d_j},& \left(19\right)\end{array}$

using FMM and then taking FFT to obtain I(t_j) where u_j and d_j have been defined just before Eq. (18). Note that u_j has been obtained in the time independent calculation. If (1/ω_T)^j and d_j in Eq. (19) are uniformly distributed on the complex plane, the FMM can be done much faster. However, as shown in FIG. 6 and FIG. 7, the unit distribution of (1/ω_T)^j and d_j are highly non-uniform in this case. Actually, (1/ω_T)^j are distributed non-uniformly along the circle (FIG. 6) while d are distributed in a sector of unit circle (FIG. 7). This makes the calculation more difficult. For N=10⁴ and T=N⁸, it was determined that the optimum number of levels in FMM is 10. With 10 levels in FMM, over 60% of the CPU time was spent on the direct sum in the FMM calculation. Using a single core, it took 3,400 s to finish the time dependent part of calculation and about 90% of the CPU time was spent on FMM. For comparison, it took 33,800 s to perform the time independent calculation. This shows that the algorithm II is of order O(1) if T<N².

FIG. 6 shows the testing of the algorithm IIa, which is suitable for T=N and the algorithm IIb designed for T=N², against the algorithm I For T<N, the results of the algorithm I, the algorithm IIa and the algorithm IIb are on top of each other. For T>N, the calculation was done for T=N². FIG. 6, only shows the results for T<40,000. There is no significant feature in the transient current plot beyond that.

FIG. 7 shows the distribution of d_j on the complex plane.

FIG. 8 shows the transient current calculated by Algorithm I, Algorithm IIa and Algorithm IIb, IIa and IIb refer to the cases with FMM methods targeting T=N and T=N², respectively.

In summary, the exact solution of the algorithm for the transient current always contains time dependent parts. According to the present invention the expression for the current is separated so that, e.g.: for a function f=f(t,x) that depends on time t and space x, it is separated into two part so that f=g(t)h(x). Thus, h(x) can be calculated first, which is a very complicated function involves many multiplications. Eventually, for any time t, the h(x) only needs to be computed once.

In addition, the major difference between Algorithm I and II is that for the expression f=g(t)h(x), in algorithm I the multiplication is performed directly; but, for algorithm II, FFT and FMM are adopted to further speed up the multiplication between g(t) and h(x).

While the present invention has been particularly shown and described with reference to preferred embodiments thereof; it will be understood by those skilled in the art that various changes in form and details may be made therein without departing from the spirit and scope of the invention.

Claims

1. A method of transient current evaluation to determine the response time and peak transient current of nano-chip in nano electronics, comprising the steps of: inputting Hamiltonians, Heq, Hneq;preparing self-energy of leads by using complex absorbing potential (CAP);generating eigenstates and eigen-energies;expanding Fermi functions from the Pade spectrum decomposition (PSD) method;constructing the exact equation of transient current which is beyond the wideband limit (WBL) from the eigenstates and eigen-energies and the Fermi function based on non-equilibrium Green's function (NEGF) formalism using residue theorem;separating the expression into space dependent and time depending components respectively so that it is of a general form I(t)=1t(Io, Vm(t), Mn), such that it is a function composed of time dependent Vandermonde matrices Vm and other space dependent matrices Mn;computing and constructing space dependent matrices Io and Mn by an optimized matrix multiplication process,compute the multiplication of Vm and Mn directly;summing up all of the contributing parts of I(t) from the multiplication of Vm and Mn and Io to form Algorithm I, which has a computational complexity of 50N3+TN2, where N is the system size and T is the time step; andoutputting the complete transient response I(t) over the user-defined time period.

2. The method of claim 1 further including the steps of: using FMM and FFT methods to compute the multiplication of Vm and Mn;summing up all of the contributing parts of I(t) from the multiplication of Vm and Mn and Io to form Algorithm II, which has a computational complexity of 50N3+2N2 log2 N for T<N2 and large N; andoutputting the complete transient response I(t) over the user-defined time period.

3. The method of claim 1 wherein the Hamiltonians, Heq, Hneq are derived from first-principles, a tight-binding method or a similar technique.

4. The method of claim 1 wherein the first principles of the input step are (DFT), TB,

5. The method of claim 1 wherein the self-energy of leads is from CAP.

6. The method of claim 2, wherein T=N and Algorithm IIa is formed.

7. The method of claim 2 wherein T<2N and Algorithm IIb is formed.