Non-equilibrium Green function method: theory and application in simulation of nanometer electronic devices

We review fundamental aspects of the non-equilibrium Green function method in the simulation of nanometer electronic devices. The method is implemented into our recently developed computer package OPEDEVS to investigate transport properties of electrons in nano-scale devices and low-dimensional materials. Concretely, we present the definition of the four real-time Green functions, the retarded, advanced, lesser and greater functions. Basic relations among these functions and their equations of motion are also presented in detail as the basis for the performance of analytical and numerical calculations. In particular, we review in detail two recursive algorithms, which are implemented in OPEDEVS to solve the Green functions defined in finite-size opened systems and in the surface layer of semi-infinite homogeneous ones. Operation of the package is then illustrated through the simulation of the transport characteristics of a typical semiconductor device structure, the resonant tunneling diodes.


Introduction
For five decades of development, the 'non-equilibrium Green function' (NEGF) formalism has become a useful and powerful calculation/computation method to the study of dynamical processes in non-equilibrium many-body systems. Historically, the NEGF theory stems from efforts in applying ideas established in the quantum field theory to treat problems in statistical physics. Accordingly, the expectation values of observables and the correlation functions between observables can be systematically determined if knowing Green functions which are appropriately defined. For instance, the correlation function between the two observables a and b measured in a system wherein a ground state Φ , 0 or a thermal-equilibrium state are established, can be calculated through the two-time Green function defined by [1][2][3][4][5]  For the progress of the field, one should mention to the key contributions of Martin and Schwinger in the formulation of the general formalism of the N-particle Green function [6]. However, in the practical point of view, in order to calculate such N-particle Green functions one needs to introduce appropriate approximations, which usually lead to the appearance of many novel correlation functions. In 1961, Kadanoff and Baym discussed some fundamental aspects in the implementation of the Martin-Schwinger theory in the view point of macroscopic conservation laws, the continuity equation for instance [7]. The development of the quantum field theory for thermal equilibrium many-body systems led to the formulation of the temperature Matsubara Green function theory. This theory is based on the concept of imaginary time, which was introduced based on the analogy between the time evolution operator and the equilibrium statistical one [2][3][4]. However, different from the equilibrium situation, non-equilibrium systems do not have neither a ground state nor a unique statistical operator. Treatment of such systems therefore requires the introduction of novel Green functions to quantify typical statistical properties such as the distribution of particles on the spectrum of microscopic states. In 1965, Keldysh introduced a way to unify a class of four basic two time-variables Green functions using the concept of timecontour [8], defined on which is the time-contour Green function or the Keldysh Green function. Inversely, by making the continuum limit, which is expressed through the so-called Langreth theorem [9], the time-contour Green function leads to a set of four real-time two-point Green functions, which simultaneously describe the two fundamental aspects of many-body systems: the spectrum of quantum states of particles and the state populations or the average number of particles occupying each state. These real-time two-point Green functions are thus able to describe generally manybody systems without needing the distinction of their statistical state, i.e., equilibrium or non-equilibrium. Since 2000, a series of conferences had been held to report progresses and to sketch perspectives of the NEGF theory. Collections of papers published in a series of proceedings named 'Progress in Non-equilibrium Green functions' report interesting key results and breakthroughs in applying this method to various fields such as plasma physics, solid state theory, semiconductor optics and transport, nanostructures, nuclear matter and even high energy physics [10].
In the last decade, the NEGF method has been developed to investigate and predict transport properties of nanoscale materials as well as the operation and the performance of nanoscale devices. The NEGF method has been implemented into many computer packages with different sophisticated levels, for instance, combined with density-functional calculations [11][12][13], or based on effective physical models [14][15][16][17]. However, many technical and physical issues still need to be solved to improve the effectiveness and the performance of such computational tools in various applications. Our computer package named OPEDEVS, standing for 'Opto and Electronic Device Simulation' [18], has been recently developed on the basis of the NEGF method to the study of the electronic and transport properties of nanoscale materials and devices with different geometrical structures. The purpose of this paper is therefore to present a short review of fundamental aspects in the implementation of the NEGF method into the OPEDEVS package. Nevertheless, all basic concepts and key equations of the theory are presented in detail, in a more suitable form for the computational purposes rather than for a full description of the theory.
The layout of the paper is presented as follows. In section 2, we review progresses in the last decades of the NEGF method in nano-electronics. In section 3 we present the definition of the four real-time Green functions, viz the retarded, advanced, lesser and greater ones, and their equations of motion and basic relations among them. In particular, we discuss the issue of how to calculate important physical quantities from the Green functions in the quantification of the transport properties of nano-electronic systems. In section 4, two efficient recursive algorithms are presented, which are implemented in the OPEDEVS package, to solve the Green functions defined in finite-size opened systems and in the surface layer of semi infinite-size homogeneous systems. Section 5 is dedicated to the illustration of the application of OPEDEVS to study a typical semiconductor device structure, the resonant tunneling diode (RTD). Section 6 concludes the paper.

Progress of NEGFs method in nano-electronics
Since the invention in 1962, the technology based on the complementary metal-oxide-semiconductor (CMOS) circuit, named the CMOS technology, has been rapidly and strongly developed in the trend of minimizing the size of the building blocks such as the p-and n-type metal-oxide-semiconductor field-effect transistors (MOSFETs). This trend, however, for the last decades, revealed severe obstacles generally known as the short-channel effects, when the device gate length approaches the nanometer scales. Though the 30 nm technology (implying that the typical length of the gate electrode of the MOSFET structure is about 30 nm) had been started since 2007, the International Technology Roadmap for Semiconductors [19] predicts that this trend will be soon stopped, because of the non-preservation of the transfer characteristics as those of conventional device structures. Fundamentally, such problems originated from the emergence of quantum effects that limit the device scaling trend.
Therefore, it has been crucial to study quantum behaviors and transport mechanisms of charge carriers which govern the characteristics of nano-scale devices. On the one hand, that research helps clarifying the lower limit of the scaling process. On another hand, it may shed light on the finding of alternative fabrication methods, and even on the invention of novel device generations. In such context, it is necessary to develop efficient theoretical approaches and powerful computational tools to study rigorously the role of quantum effects involving in the transport of charge carriers in nanoscale device structures, and hence their dominant influence on the device performance. Figure 1(a) illustrates a typical device structure of two terminals. The central region is called the active one that governs the behaviors of the whole device, while the two end regions, called the leads, play the role of the source and the sink of charge carriers. The active region is thus considered as an opened system which communicates with its environment via the leads. At the scales of conventional device structures, with the gate length in the order of μ 1 m, various classical and semi-classical approaches, for instance, the drift-diffusion model and/or the ones based on the semi-classical Boltzmann equation, have been developed [20]. Such approaches, however, become invalid in capturing the quantum effects, such as the quantum interference, confinement and tunneling, which become emerging when the size of the devices is comparable to, or even smaller than, the phase breaking length ϕ L of charge carriers. When the wave function of carriers can extend largely over the central region to the leads, carriers can move coherently through the whole device. The transport, in this sense, is said to be in the coherent regime, or in the ballistic transport regime. The current of charge carriers flowing through the structure therefore can be determined by a rather simple and intuitive formula introduced by Landauer in 1947, which was then validated and generalized for multiprobe systems by Buttiker using the framework of the scattering theory [21][22][23][24], where e 0 is the absolute value of the electron charge, and  ( ) E is called the transmission coefficient, which is actually the summation of transmission probabilities (the transmittances) of charge carriers over all possible single-particle modes. The transmission coefficient is therefore a physical quantity quantifying the 'transparency' of the opened system, which is viewed as the agent scattering the carrier plane waves incident from the terminal regions. In the above equation, the role of the carrier source (the left-L) and sink (the right-R) regions is explicitly reflected through the presence of the two functions ( ) n E L and ( ) n E R which define the average occupation number of carriers on states with energy E. Since these regions are usually much larger than the device active one, they are usually assumed to be in a thermal equilibrium state, and thus characterized by a temperature ν T and a chemical potential are usually approximated by the Fermi--Dirac function with some special notices involving in the specific value of the temperature and the chemical potential [25]. The calculation of the current of charge carriers by equation (2) therefore revolves around the determination of the transmission coefficient. Different methods may be used to this aim, including the use of the scattering method with the help of the transfer matrix technique [26], or the use of the real-space Green function method associated with the Fisher-Lee relation [27].
Though equation (2) is useful, it is based on the singleparticle picture. In practice, even when the ballistic and coherent transport regime dominate in the nanoscale structures, the problem of carriers scattering off impurities, lattice vibrations, and the surfaces/interfaces in/of material layers is in general inevitable. These scattering processes may even govern typical behaviors of devices, for instance, the effect of phonon-assistance transport and the super-Poissonian noise in double-barrier resonant tunneling structures [28]. Such effects apparently cannot be microscopically accounted for in the framework of the Landauer-Buttiker formalism (the scattering theory). A general theory that enables naturally connecting the two transport regimes, i.e., the ballistic and the diffusion ones, is therefore needed.
To the simulation of electronic device structures, the transport equations of charge carriers, in principle, must be presented and solved in the real time-space. In 1990, Caroli pointed out that some physical quantities quantifying the transport characteristics of an opened system can be expressed in terms of the Green functions which are defined only in the active region [29]. In 1992, Meir and Wingreen based on the NEGF formalism to formulate an expression for the charge current flowing through a finite scattering region and recovered the well-known Landauer formula in the limit of no-scattering, i.e., in the coherent transport regime [30]. In the next sections, we will review these key results since they play the basic role in the development of the NEGF formalism for the device simulation. In 1995, Datta introduced a book wherein many practical aspects of implementing the NEGF method are discussed in a rather pedagogical way using the matrix representation of the Green functions of the two terminal structures [31].
Practically, the most difficult part in the use of the NEGF method for investigating transport properties of a system lies in the treatment of many-particle interaction effects. Diagram techniques developed for the many-body Green functions in general can be used in the NEGF formalism. It, however, is a formidable task [4]. For the cases in which interaction processes play the role of making the transition of particle from one state to another, but do not cause significant changes in the dispersion relation, the NEGF formalism can be used in the quantum-version form of the Boltzmann equation for the so-called Wigner function [32]. Mathematically, the Wigner function plays the role similar to the Boltzmann distribution function in the determination of the kinetic properties of a quantum system, but it is not positively defined. In addition, since the transition processes can be treated using the framework of the Fermi golden rule, many useful methods and techniques already developed in the classical statistics, such as the Monte-Carlo method, can be used to solve the quantum Boltzmann equation [33,34].
Regarding the solution of the NEGFs, it is worth noting their recursive structure in terms of the so-called Dyson equation. The fact is that the realization of a mathematical structure in a system of equations will strongly facilitates computational processes. This point will be detailed in section 4, but informatively, it should mention some important contributions in the field of numerical calculations. In 1997, Lake et al. successfully used the recursive method to calculate all diagonal elements of the retarded Green function matrix defined in the device active region [14]. This technique enables combining different parts of a whole device using the concept of coupling self-energy. In 2002, Svizhenko et al. started from the matrix structure of the transport equations to introduce a rather general algorithm to compute not only important elements of the retarded/advanced Green functions but also of the other correlation functions, i.e., the lesser and greater Green functions (see next section) [15]. This algorithm is very useful and now commonly used to find, in a systematic way, all of the Green functions defined in the device active region given the self-energies characterizing the terminal (source/sink) regions. For some simple transport models, such coupling self-energies can be analytically found [26,31]. However, when dealing with complicated models, such as multi-bands and/or atomic tight-binding description, it is not easy to handle the analytical calculations. The algorithm proposed by Sancho and Rubio in 1984 to treat semiinfinite systems was usually revisited. It is now known as a rather general, stable and efficient method to calculate the coupling self-energies [35]. All these techniques are the key ingredients for the implementation of the NEGF method to the problem of quantum device simulation. In the next sections we will review these points, which have been successfully implemented in the OPEDEVS package.

Definitions
Given an electron system which is described by a Hamilto-nianĤ represented in terms of the field operatorsˆα β † c c , , where α and β denote, for short, the sets of quantum indices labeling quantum states of the system. The time-contour Keldysh Green function is defined by (see appendix A): , wherein  denotes the time-contour such that  τ τ ∈ , 1 2 (see figures 2(a) and (b)),  T is the time-ordering operator on the time-contour, and the angle bracket denotes the ensemble average as specified in equation (1). The time-contour  can be seen as composed of three branches The definition of the NEGF given by equation (3) actually is the general one for the four real-time Green functions: the lesser < ( ), greater > ( ), causal (time-ordered) c ( ) and anti-causal (anti time-ordered)c ( ) Green functions when taking the limit η → 0 [2][3][4]: where the curly bracket denotes the anti-commutator. From equations (5)-(8) we clearly realize the general relationships among the four real-time functions: In most cases, we are usually interested in the steady state of a system. The two-time correlation/Green functions become dependent on only the difference between the two time variables¯= − t t t . 1 2 Therefore, it is more convenient to work with their Fourier transforms α β Inversely, we have 1 2 In the above equations and so later, the superscript ν is used to denote the ones ≶ { } r a , , (except some cases which will be noted), and E is the energy variable. Practically, it is useful to remark the dimension of the defined Green functions in equations (4) and (13). Since the field operators are usually chosen to be dimensionless, the real-time Green functions and their Fourier transforms α β in 'the International System of units', therefore have the dimension of (J s) −1 and J −1 , respectively.
For electron systems, it is worth to mention to the spectral density matrix α β ( ) which is defined as follows , , , , As generally presented in [2,3] the diagonal elements, i.e., when α β = , of this quantity are positively definite and provide information of the energy spectrum of elementary excitations in the systems. The spectral density matrix therefore has the dimension of the energy inversion.

Equations of motion
In appendix B we present the Dyson equations for the timecontour Green function  τ τ′ ( ) G , (In the following the bold characters are used to denote the matrix and vector notation, i.e., to omit quantum indices). Mathematically, it is an integral equation obtained from the perturbation expansion of equation (A.12). In this section, we will apply the Langreth theorem (proved in appendix B) to the Dyson equations (B.5) and (B.6), to specify the equations of motion, in the general form, for the four real-time Green functions, i.e., the retarded, advanced, lesser and greater.
Indeed, making use of the rules given in table B1, the continuum limit of equation (B.5) becomes [32]: r a r a r a Similarly, the continuum limit of equation (B.6) yields r a r a r a Because of the time-convolution integrals in the above equations, it is useful to write down the equations for the Fourier transforms of the Green functions. Straightforwardly, it results in r a r a r a 1 like to emphasize that all the above equations for the four realtime Green functions are exact and general in the meaning of the perturbation expansion, i.e., when one can define the self-energies. Given an expression for the full Hamiltonian =ˆ+Ĥ H H, 0 1 the headache problem always is how to specify the self-energies. Technically, the self-energies can be constructed using the Feynman diagram technique [2,3], or the method of equation of motion. These two are the most general, common and powerful techniques. For the latter, readers could consult [36] as an example.

Some basic properties and relations
In this section some general analytical properties of the Green functions are summarized. From the definition of the real-time Green functions, it is straightforward to check that where E(k) expresses the dispersion law of electron, and V the volume of the system. Substitute these expressions into the definition of the lesser Green function, we obtain is the average number of particles occupying the state with energy ( ) . The Fourier transform of this function is thus easily determined to be: In the same way, the expressions for the greater, retarded, and advanced Green functions are derived as follows r a Next, from equation (17) we can show that the retarded/ advanced Green functions are analytic in the upper/ lower one-half complex plane, respectively (see equation (36) as an example). This property allows expressing G r/a (E) in the form: r a where η is a positive infinitesimal number. Equation (37) is well known as the Kramer-Kronig relation of analytic functions [2,3]. In general, the self-energies are complex functionals of the Green functions. However, in some cases, for instance, the coupling of an open system with its leads (see equation (57)), and/or the phonon-electron interaction in the self-consistent Born approximation [14], the self-energies simply take the linear form as follows could be the device-lead or electron-phonon coupling factors. In these cases, we can write down an equation similar to equation (37) for the retarded/advanced self-energies: Similar to the spectral density matrix ( ) A E , the diagonal elements of Γ ( ) E also have a physical meaning. Specifically, they are determined as the one-half widths of spectral peaks [3].
In opened systems, we distinguish two kinds of sources which induce the changes in the spectral and kinetic/transition pictures of charge carriers: one involves in the connection of the opened system to its environment (the device-lead coupling), and the other the intrinsic interaction/scattering processes of charge carriers inside the system. The self-energies (17) It can be shown (see section 3.5) that when assuming the leads as semi-infinite homogeneous systems, the lesser and greater self-energies Σ ≶ c explicitly take the following form [3,32]:   (40) and (57).

Electron and hole densities
The aim of this section, and of the next one, is to make a link between the Green functions defined in the previous sections and some physical quantities which play an important role in the transport problem. Such physical quantities include the distribution densities of charger carriers and the charge current density. We will show that knowing the Green functions is sufficient to calculate these physical quantities.
In order to formulate in detail all the following expressions in our transport problems we need to specify the quantum indices α and β labeling the operators and Green functions. In our convention, these quantum indices refer to the sets of the three ones λ { } q k , , , in which the first one, denoted by Latin characters such as p or q, is used to label the unit cells or the meshing nodes along the transport direction, see figure 1(b). The second index λ labels the conduction electronic bands or the atomic valence orbitals contributing to the transport in a unit cell. The last one k is an option, which may be used depending on the specific geometrical structure of the systems under consideration to characterize the transverse motion of particles. When the cross-section of a device structure is sufficiently large so that the effect of its boundaries is not important, k serves as a wave-number vector. The illustration of using these three indices is presented in the sketched device structure in figure 1(b).
Denote η s the factor counting any degeneracy of some quantum states (spin and/or valleys of energy bands, for instance) the electron density distributed in a unit cell, or in a mesh-spacing, is determined as the expectation value of the electron number operatorˆ=ˆλ where in S is the cross-section of the system and a x is the width of the unit cell or the mesh-spacing. Using the definition of the lesser Green function given in section 3.1 we can express n q e as follows: By the same manner, but noting to the anti-commutation relation between the creation and annihilation operatorŝˆ= we can obtain a similar expression for the hole density in terms of the greater Green function, Knowing n q e and n q h the distribution of charge density in a system is determined, A are the donor and acceptor densities as the ground charges. The charge density is necessary to describe the electrostatic effects which affect on the formation of electric current in the system. This point will be clarified in section 5.

Charge current density
We show that the charge current density can also be calculated by means of the Green functions. To do so, it is necessary to find an appropriate expression for the electric current density in terms of the basic operators. Because the formation of a charge current is the result of the response of a system to a transverse electromagnetic field, which is characterized by a vector potential A(r) [37]. Remember that the Hamiltonian of a system, in the second quantization representation, can be written as the summation of a singleparticle term and other many-particle ones [4,5], for instance, wherein the 'hopping' integrals H p q , is a functional of the vector potential according to the Peierls substitution, and is given by (see appendix C): pq p q The velocity, and thus current, operator is thus determined as follows [4,38]: In the approximation of nearest-neighbor coupling this formula becomes x q q q q 1 1 By making the ensemble averaging of this operator, with the notice of the definition of the Green functions, we obtain the expression for the density of electrical current flowing through the cell q as follows [14,15]: This formula is general and one can use it to calculate directly the current density. However, since the current density along the transport direction is uniform in the steady state, i.e., = = I I I, p q it is interesting to reformulate equation (52) in a rather convenient form. Indeed, by considering the cell = q 1, the matrix blocks < G 10 and < G 01 of the lesser Green function matrix are specified as follows: It is helpful when noting to the fact that the leads contacting to the device active region are usually large, compared to the active region. They therefore could be assumed to be in a (quasi-) thermal equilibrium state, i.e., being characterized by a function η n x ( ) determining the average occupation number of charge carriers on a quantum state (η = L R , labels the left and right lead regions). According to equations (34)-(36) the lesser Green function matrix cL r a r a 1,0 0,0 0,1 the retarded/advanced (left) lead-device coupling self-energies and use equation (40), it yields Substitute this result into equation (52), the current density is thus rewritten in the form This formula for the current density is exactly the result derived by Meir and Wingreen in reference [30]. We can go further by noticing that  (41)). Substitute these expressions into equation (58) we are able to separate the current density I into the two contributions I c and I , s in which    (64) is therefore interpreted as the probability of particle, occupying the transverse mode k, transmitting through the device active region. In the case of equilibrium states, both current density components I c and I s become vanished automatically since For the latter component, i.e., I , s it is clearly induced by the unbalance of interaction/scattering processes of charge carriers. This term clearly does not contribute to the total current density if interaction/scattering effects are not taken into account.

Numerical methods and techniques
In the previous section, the key concepts and equations of the NEGF method have been presented as the basis for employing the theory as a calculation method. Given a Hamiltonian and functional expressions for the self-energies, it looks straightforward to just perform some matrix algebras, such as the matrix inversion and the matrix multiplication, to obtain the Green function matrices, see equations (19) and (20). It is, however, in practice, a formidable task because all involved matrices are usually very larger. Therefore, the numerical computation is very expensive in the meaning of consuming huge computer resources. Fortunately, as shown in the previous section, only a small number of elements of such matrices is needed to calculate physical quantities, for instance, the charge carrier densities and the charge current density, see equations (44), (45), and (52) or (62) and (63). It is therefore very useful to develop numerical methods such that only interested matrix elements are calculated. In this section, we will present two such algorithms, one for solving the Green functions defined in the finite active region, and the other for the Green functions defined only on the surface of semi-infinite homogeneous regions.

Recursive algorithm for the retarded/advanced Green functions
We first present the algorithm for solving equation (21), i.e., calculating some important elements for the retarded/ advanced Green function matrices. Since the equations for the retarded and advanced Green functions have the same structure, the following presentation is for the retarded one.
Given a Hamiltonian H 0 and a self-energy Σ , r which are represented as tri-diagonal and diagonal block matrices, respectively, to characterize the dynamics of charge carriers in the device region. Call N cell the number of matrix rows/columns. From equation (21), the retarded Green function matrix is determined as the inversion of the matrix To invert the right-hand side matrix we separate it as follows The recursive structure of this equation suggests the following view: the whole device active region can be seen as the construction of N cell cells stacked sequentially, starting from the first cell = q 1 (left-contacted), or from the last one = q N cell (right-contacted). Accordingly, let us consider a system of +  Equations (73) and (74), with the notice of equation (75), constitute the so-called 'forward' calculation procedure for solving the retarded Green function of the opened system. The 'backward' procedure is therefore necessary to complete this task. To do so, we should specify the meaning of the 'bared' function G r0 . Assume that we already know the block + + G q q r 1, 1 . The retarded Green function of a system of (q + 1) cells, but the (q + 1)th cell does not couple to the block of q Combine these, we obtain the following recursive equation: The 'forward' and 'backward' calculation procedures as described can be directly applied to calculate the advanced Green function. However, in practice, this function is determined by invoking the relation ⎡ ⎣ ⎤ ⎦ = † G G a r rather than performing again the recursive algorithm.

Recursive algorithm for the lesser/greater Green functions
A recursive algorithm can be also developed to calculate the lesser/greater Green functions as shown by Svizhenko et al in [15]. In this section, starting from the Dyson equation for these functions we will re-establish that algorithm with two typical steps of 'forward' and 'backward'. Remember that the Dyson equation for the lesser/greater Green functions is not the same as that for the retarded/advanced Green functions. In appendix B we present the rules for making the continuum limit of the time-contour Dyson equation to derive the realtime Dyson equations for the lesser/greater Green functions [2,9,17,32]. Since the equations for the lesser and greater functions have the same mathematical structure, following we only present the algorithm for the former one. From equation (22) the lesser Green function matrix < ( ) In order to find < G in a recursive manner we also construct the left-contacted lesser Green functions the left-contacted lesser Green function of the system with only q cells. Equation (86) obviously has the role similar to that of equation (73) in the construction of the 'forward' calculation procedure with the same notice that the cells with = q 1 and = q N cell couple to the left and right leads (it means that the contact-induced lesser self-energies should be appropriately integrated into the general self-energy Σ + + We further use equation (67) to expand G a , and thus obtain Accordingly, we derive the following recursive equations: These results obviously constitute the 'backward' procedure to completely determine all diagonal and some important off-diagonal blocks of the lesser Green function matrix. For the greater Green function, since there is no simple algebra relation to the lesser and/or retarded/advanced functions, it must be independently determined. Fortunately, since both ≶ G have the same mathematical structure, the algorithm described by equations (86)-(92) can be directly used to calculate > G .

Sancho-Rubio algorithm
As shown in the previous sections, to find the Green functions of an opened system one has to already know the self-energies characterizing the coupling between this system and its environment. Given H 0,1 and H 1.0 the Hamiltonian matrix block defining the coupling between the cell = q 1 of the device active region and the surface cell of the left lead. According to equations (57), (41) and (42) we see that it is mandatory to determine the block g rL 0,0 of the left-lead retarded Green function. In the case that the lead regions are finite in size, the recursive algorithm presented in section 4.1 can be directly used. However, because the lead regions are usually much larger than the active region, they could be seen as infinite-size systems. This assumption actually can simplify the calculation for g rL 0,0 . Indeed, the block g rL 0,0 can be formally written in the form Practically, iteration schemes can be used to solve directly equation (94). However, the convergence is usually very slow. In 1984, Sancho and Rubio pointed out that the convergent rate can be strongly improved by taking special care of the recursive structure of the Dyson equation. The scheme presented below is for such care and now known as the Rancho-Rubio algorithm [35].

M E is called the transfer matrix and is recursively given by
Denote g rL the left-contacted retarded Green function of a semi-infinite system. By invoking Dyson equation (67) we obtain the following results for the blocks g , rL r 0,0 0,0 1 0,1 1,0 and for ⩾ q 1,   Accordingly, we have The process is repeated until + t n 1 and ′ + t n 1 as small as needed, then ≈ + g 0

OPEDEVS package and the transport characteristics of RTDs
In this section, we present some calculation results obtained from a numerical investigation of the transport characteristics of a typical electronic device, the RTDs. The calculation was performed using module HETS integrated in the OPEDEVS package [18], whose construction is based on the NEGF method described in the previous sections. The module HETS (appreciation of 'Hetero-semiconductor structures simulation') was designed to solve the electronic transport problem of hetero-semiconductor structures. Specifically, HETS is a collection of computer programs to numerically solve a couple of two differential equations, the effective Schrodinger equation and the Poisson one, which are defined in a finitesize active region of length L and a sufficiently large crosssection. These two equations respectively read is the principle schema of the RTD. Accordingly, this device can be seen as made up of five principle semiconductor layers. The center layer plays the role of a quantum well (W). It is sandwiched between two other semiconductor layers as the potential barriers (B1, B2). The system of these three material layers is then connected to two other thick semiconductor blocks, which are doped to make two good contact regions, the emitter (E) and the collector (C). The doping is not necessary uniform. In practice, the dopant concentration in the outside-end regions of these two layers is high, while it is very low in the regions close to the center one. Regarding this fact, one can distinguish the buffer regions (BF1, BF2) from these contact regions. In figure 3(b), we draw the profile of the conduction band edge in the flat form to illustrate the formation of the potential barriers due to the finite offset energy. As shown below, the center layer and the two adjacent ones essentially govern the transport characteristics of the whole structure. It is thus called the device active region.
The operation of the HETS module is governed by the parameters defined in three input files. The first and second files, named GEOMETRY_STRUCTURE and MATERI-AL_STRUCTURE, are designed to specify the values of the parameters which define the geometry and the materials in use of the device. The third input file, named OPER-ATION_PARAMETERS, specifies the initial values for some physical quantities as well as other parameters to control the numerical calculation, the accuracy of the self-consistent electrostatic potential, for instance.
Once the three input files are established (actually the templates are provided in the package), the first step of the calculation procedure is to perform the self-consistent solution of the Schrodinger and Poisson equations to obtain the correct values of the electrostatic potential and the carrier density distributions. This step is realized with an initial guess for the electrostatic potential. It can be made manually or by choosing the default option. The self-consistent calculation is performed by invoking the program ModelCalculation as a line-command as follows In the first form, the option '-s' is specified to confirm the string 'suffix' as the suffix of output files that the program generates to save the values of the electrostatic potential and of the carrier densities. For instance, if 'suffix' is specified by 'Vp0.5', the output-files names will be 'Pots_Vp0.5' and 'Dens_Vp0.5'. In the case that the option '-s' is used but 'suffix' is forgotten to provide, the saving mode of the program will set the output-file names in the default form as 'Pots_' and 'Dens_'. Running the program by this command, it requires to specify 'input_file' as the file name, including the path of course, containing an appropriate value for the electrostatic potential U(x) as the initial guess. This running mode can be used to restart the calculation from a previous step in the cases that the self-consistent process is stopped (by electric cutting, for instance) before reaching the convergence. The second form of the command is convenient to start the calculation process from the beginning, with the option '-d' added to imply that the default value of U(x) is used as the initial guess.
Once the self-consistent solution for the electrostatic potential and charge carrier densities is achieved, we can perform the calculation for interested physical quantities. For instance, the charge current density at a given voltage can be calculated by invoking the program 'ChargeCurrent' via the line command: where 'input_file' is a string specifying the file (including the path, of course) saving the electrostatic potential value (i.e., the 'Data/Pots_suffix' file) and 'output_file' specifies the target for saving the current value as output. In the thermal equilibrium state, the module HETS also provides a program, named 'Conductance', to calculate the electrical conductance at two regimes, the absolute zero temperature and the finite temperature ones, using the command: By varying the values of some parameters in the input files, one can change the operation condition of the simulated device. One thus can make an investigation of some transport properties of the device, for instance, the current-voltage (I-V) characteristics or the conductance. Note that when invoking the program 'ChargeCurrent', it results in only one value for the current density corresponding to the value of the bias voltage specified in the file 'OPER-ATION_PARAMETERS'. For the program 'Conductance', it also results in one value for the conductance for a provided temperature, except for the case of zero temperature, a file saving the conductance value as a function of Fermi energy. In OPEDEVS we, however, commonly use the shell-script technique to automate some calculations, for instance, automatically calculating all desired points of the I-V curve.
Apart from the described functionalities, module HETS provides some other utilities to calculate various fundamental physical quantities, which may help to analyze the physical picture behind the transport characteristics. For instance, if invoking these three commands:

SpectralFunctions
one can obtain, respectively, data for the contact-surface Green function, or the transmission coefficient, or the local density of states (LDOS) and the carrier distribution functions in the real-energy space. Module HETS is an open package in the meaning that users can freely develop different programs interfacing with the program 'ModelCalculation' to calculate interested physical quantities.
We now present the simulation results for a concrete RTD. In figure 4, we display the I-V characteristic of a RTD sample with all parameters given in table 1. The figure obviously shows the well-known N-shape of the I-V curve [16,39]. In the small bias voltage range, the current increases when raising up the value of bias voltage. However, when the voltage passes one critical value, called V c (which determines the current peak I P ), the current suddenly falls to a much lower value (the current valley I V ), and then slowly increases if continuing increasing the voltage. This typical non-linear form of the I-V curve, characterized by a very large negative differential resistance, had been a topic attracted an intensive consideration of both fundamental and technological researches in the 2000's period [14,16,40].
To shed the light on the operation of the RTDs, we present in figures 5 and 6 the profiles of the conduction electron density and the conduction-band bottom for several values of the bias voltage. In the case of no bias voltage ( = V 0), the device is in the thermal equilibrium state. Since the structure of the considered device sample is symmetrical, the profiles of both quantities are symmetrical through the center region. Comparing the profile of the conduction band bottom, which was self-consistent calculated (figure 6), to that of the energy offset plotted in figure 3(b), we obviously see the effect of the non-uniform doping in the emitter and collector as the rise of the potential in the buffer and center regions. When increasing the bias voltage, we see the down-shift of the conduction band bottom in the collector region (right). The profile of the electron density is thus changed into the asymmetrical form. Initially, when the voltage increases, we see the increase of the electron density in the quantum well. However, when the voltage passes    To find the answer, we need to compute and analyze other physical quantities. In figure 7, we show the LDOS of the simulated device for = V 0. The figure clearly shows the existence of LDOS peaks in the quantum well region. These peaks obviously reflects the quantization of the electron states in the quantum well due to the confinement caused by two finite potential barriers. We found that the position of such quantized levels, relative to the conduction band bottom, are very consistent with the values estimated from the equation for finite barrier quantum wells [39], According to the Landauer-Buttiker formula (equation (2)), the transmission coefficient is a quantity characterizing the transport of electrons. In figure 9, we present the results for this quantity obtained by invoking the program 'CheckTransmission'. The figure was plotted for three curves corresponding to = V 0.00, = = V V 0.26 c , and = V 0.34 V. The transmission coefficient for each value of the transverse vector k is nothing rather than the transmittance of electron through the system. The figure, plotted for = k 0, apparently shows a sharp peak in the energy range lower than the potential barriers, and an oscillation behavior in the energy range above. The peak was found to point exactly the position of the first quantized level in the quantum well. It therefore implies that the electron in the emitter region can easily tunnel through the double potential barriers when its energy aligns with the quantized level in the quantum well. When increasing the voltage, this transmission peak shifts down due to the band bending and then disappears when V surpasses V c . The disappearance of the transmittance peak is simply due to the position of the quantized level in the   quantum well does not match with any conduction states in the emitter region. The quantized level thus becomes lost its role of a coherent conducting channel.
To analyze the formation of the charge current, we additionally present in figure 10 the spectrum of the current density I(E). This quantity is defined as the function of energy under the integral in equation (2). The calculation was performed by invoking the utility 'CurrentSpectrum'. As expected, the current spectrum shows the significant peaks corresponding to the transmittance peaks with different transverse modes k. These behaviors apparently are the evidences confirming that the transport of electrons in the RTD device is governed by the quantization of electron states in the quantum well. An electron in the emitter region can be transmitted through the double-barrier system if it occupies a state whose energy aligns with the quantized level in the quantum well. By this way, the transport of electrons is in the resonant regime. That is why the term 'resonant tunneling' is commonly added to the device name.

Conclusion
We have presented a short review of fundamental aspects in the implementation of the NEGF method into the computer package OPEDEVS that we have recently developed to study transport properties of nano-scale device and low-dimensional material structures. The definition of the four real-time Green functions, the retarded, advanced, lesser and greater Green functions as well as their equations of motion and basic relations among them are presented in the matrix form, which is suitable for the numerical computational purpose. Particularly, by rewriting the equations of motion of the Green functions in the form of the Dyson equations we have reconstructed two efficient computation schemes, one for solving the four Green functions defined in finite-size opened systems and the other for finding the Green functions defined only in the surface layer of semi-infinite homogeneous systems. These two advanced computational techniques have been successfully implemented in the package OPEDEVS. As an illustration for the operation of OPEDEVS, we have presented the computational procedure of a numerical investigation of the electronic transport characteristics of a typical semiconductor device structure, the RTD. The RTD simulation was realized by the module HETS integrated in the package OPEDEVS, which was generally designed to solve the transport problem of semiconductor hetero-structures. The transport characteristics of the device, including the I-V curve, the charge carrier density distribution, and other physical quantities, have been shown and analyzed to clarify the physical picture behind the operation of this device structure. () to make a link between the state Φ t ( ) I to the non-perturbation one Φ 0 and another one Φ + ∞ ( ) as follows We now should distinguish the difference between the equilibrium and non-equilibrium problems. For the former one, the system will asymptotically return to the initial unperturbed state as t increasing, i.e., Φ Φ +∞ ∝ ( ) 0 . This fact is well-known as the Gell-Man-Low theorem. It, however, does no longer hold for the situation of non-equilibrium. Despite of that, we still establish the relationship between Φ 0 and Φ +∞ ( ) thanks to the use of the Smatrix, , making use of equations (11)-(12) one can realize that there are more than one expression for these functions, of course they must give identical final results. For instance, if using the relation , , r c 2 2 2 , we will have:

I t t I t t I t t A t i t i B t i t i I t t A t t B t t A t t B t t A A B B A B A t t B t t A t t B t t A t t B t t
, , ,