Effective Simulations of Electronic Transport in 2D Structures Based on Semiconductor Superlattice Infinite Model

Effective simulations of semiconductor superlattices are presented in the paper. The simulations have been based on the Wannier function method approach where a new algorithm, inspired by Büttiker probes, has been incorporated into determining the Green function procedure. The program is of a modular structure, and its modules can either work independently, or interact with each other following a predefined algorithm. Such structuring not only accelerates simulations and makes the transport parameters possible to initially assess, but also enables accurate analysis of quantum phenomena occurring in semiconductor superlattices. In this paper, the capabilities of type I superlattice simulator, developed earlier, are presented, with particular emphasis on the new block where the Fermi levels are determined by applying Büttiker probes. The algorithms and methods used in the program are briefly described in the further chapters of our work, where we also provide graphics illustrating the results obtained for the simulated structures known from the literature.


Introduction
Semiconductor superlattices (SL) are modular structures where nanometre layers of two different semiconductors or their alloys are arranged alternately. They generate a periodic electrostatic potential, of which period (d) significantly exceeds the period (a) of the potential originating from atoms within the crystal lattice (see Figure 1).
Electrons moving in a SL are subject to quantum effects, so they can be used to construct innovative electronic devices, such as quantum cascade lasers (QCL) [1][2][3][4]. Despite the simple spatial geometry of the SL, the electron transport modelling in such structures has proven to be very difficult, as it requires considering many factors, such as the influence of the potentials from superlattice and the applied electric field alike, the impact of potentials derived from dopants, as well as interfaces' roughness, and last but not least, electron scattering on both photons and phonons.
In the primary models, a constant effective mass of electrons within the whole superlattice field [5,6] was assumed and its value was set as the average value of all effective masses in the structure individual layers. The potential of the SL was consistent with the Kronig-Penney model and it took the form of a rectangular wave, while the method of bonding wave functions and their derivatives at the boundaries of superlattice layers were used to determine the electron states. Over the years, the models underwent many improvements, of which the approach to solve the Schrödinger equation with a rate equation method (RE) [7,8] was particularly interesting. This method assumes the boundary conditions to have taken the form of hard walls at the edges of the structure composed of three consecutive periods of SL, and the solutions obtained for the central period to be representative of the whole method assumes the boundary conditions to have taken the form of hard walls at the edges of the structure composed of three consecutive periods of SL, and the solutions obtained for the central period to be representative of the whole system. An extension of the RE approach is the Monte Carlo method (MC) [9,10] based on solving the Boltzmann equation. With this method, it is possible to determine the system parameters both for the stationary state, and for its evolution over time. A comprehensive approach to model superlattices was proposed with non-equilibrium Green's function formalism (NEGF) [11]. Self-consistent solving Dyson, Keldysh and self-energy equations in the real space [12] (RS) can provide full information about the system, though it is computationally demanding [13][14][15]. Thus, another effective approach has been developed, namely the Wannier function method (WFM), where computations in the complex space and Wannier function base are required [16][17][18][19]. Here, the Hamiltonian matrix is small in size and its elements depend solely on the energy. This paper presents effective simulations of semiconductor superlattices, based on the WFM approach with Büttiker probes applied. This approach has neither been used to date, nor has it been described in the literature. To serve this aim, a computer program (superlattice simulator-SLSim) has been developed, where the Hamiltonian device matrix is transformed from an energy-position representation to a pure energy representation. It has resulted in a significant reduction in the matrix size. The implemented NEGF formalism allows accurate calculations with regard to basic quantum effects. The SLSim is constructed of several modules that are interoperable according to a predefined algorithm, but they can also work independently. Due to such 'structuring', initial assessment of the transport parameters for the given voltage are quick. Moreover, by combining properly selected program blocks, accurate calculations with the set of electron scattering taken into account are also possible. The numerical methods applied in each program block have been described successively in our previous works [20][21][22]. A successful comparison of this program with the exact RS model and measurement results was also carried out [23], and it confirmed our methods to be correct. This allowed us to subject the SLSim to final tests, whose results are described here.
(a) (b) Figure 1. A part of a superlattice structure made as the system of two semiconductor layers: (a) the periodic potential of d period that arises in the material as a result of semiconductor layers modulation arrangement; (b) the periodic potential derived from the crystal lattice atoms of the superlattice.
The structure and algorithms of the performed SLSIM are described in further sections of this work. All sections and methods are illustrated with simulation results to present how the program is to be used for specific superlattice structures described in the literature. We start by showing the main algorithm of the program and a list of the investigated superlattice structures. Then, we describe the method for solving the Schrödinger equation and constructing quantum states that represent three superlattice periods. The NEGF algorithm with a detailed description of the Fermi levels determination method by applying Bütiker probes is given in the next chapter. This chapter Figure 1. A part of a superlattice structure made as the system of two semiconductor layers: (a) the periodic potential of d period that arises in the material as a result of semiconductor layers modulation arrangement; (b) the periodic potential derived from the crystal lattice atoms of the superlattice.
The structure and algorithms of the performed SLSIM are described in further sections of this work. All sections and methods are illustrated with simulation results to present how the program is to be used for specific superlattice structures described in the literature. We start by showing the main algorithm of the program and a list of the investigated superlattice structures. Then, we describe the method for solving the Schrödinger equation and constructing quantum states that represent three superlattice periods. The NEGF algorithm with a detailed description of the Fermi levels determination method by applying Bütiker probes is given in the next chapter. This chapter also includes the simulation results for the charge and current density distributions, with Lo-phonon scattering taken into account. The effect of doping on the charge distribution is described in the chapter dedicated to solving the Poisson equation. The work also contains a short presentation of the program graphical interface and a short summary, where further directions for the program development are set out.

The Main Algorithm of Calculations
The procedure for determining the allowed minibands is begun by modeling the superlattice structure under assumed conditions. In this block, the Schrödinger equation is solved. It yields the Bloch's functions calculated for the energies grouped within the allowed minibands. The next stage is to construct the base of quantum states representing three periods of the considered structure. Depending on the polarization conditions of the device, the calculations are further carried to determine either maximally localized Wannier functions (MLWF), or Wannier-Stark functions (WSF) for the thermodynamic equilibrium and the voltage-polarized SL, respectively.
Having determined the base of quantum states, it is possible to calculate the matrix elements of the Hamiltonian device in the energy representation. The transport parameters are computed further within the NEGF formalism. The SLSim provides several options for calculating the Green functions (GF). Under a much less accurate, though very fast, approach, the constant values of electron scattering times are assumed and Büttiker probes are applied in the process of finding the self-consistent charge density. It allows the SLSim to calculate the transport parameters for a given voltage range. Current-voltage characteristics can also be obtained quickly, which makes the simulator an effective tool for initial testing semiconductor superlattices. Nevertheless, by implementing the exact procedure for calculating GF, although it significantly increases simulation time, the specific electron scattering can be included.
Doping in the selected semiconductor layers is taken into consideration by solving the Poisson equation. Figure 2 presents the main algorithm for SL simulations.
Based on the above algorithm, a series of simulations of structures known from the literature were made, the parameters of which are presented in Table 1. Table 1. Basic parameters of superlattice structures simulated at work. Structure no. 1 is described in [17], structure no. 2 can be found in the article [24] and structure no. 3 is presented in [3]. The capabilities of SLSim are presented in the following subchapters.

The Schrödinger Equation
The standard transfer matrix formulation method (TMF) [17] is used to solve the Schrödinger equation where Ψ (z) are envelope functions, V(z) is a spatially dependent total potential and m e (z) is the effective mass constant within each semiconductor layer. Solutions of Equation (1) occur in the form Electronics 2020, 9, 1845 4 of 16 where z j is the position of the structure region j, and k j (E) = 2m j E − V j /h 2 , with m j and V j as the mass and the potential in the region j, respectively. From the base functions (2), the Bloch states ϕ ν q (z) are computed, where ν is the index of the miniband and q is the Bloch vector within the range − π /d , π /d . The capabilities of SLSim are presented in the following subchapters.

The Schrödinger Equation
The standard transfer matrix formulation method (TMF) [17] is used to solve the Schrödinger equation where Ψ(z) are envelope functions, V(z) is a spatially dependent total potential and m e (z) is the effective mass constant within each semiconductor layer. Solutions of Equation (1) occur in the form where zj is the position of the structure region j, and , with mj and Vj as the The term d stands for the SL period length. A detailed description of the implementation of the method presented here was published earlier in [20]. The exemplary results for this part of simulation, dedicated to the structure described in [3], are illustrated in Figure 3 where the module (in red) and the phase (in blue) of the selected Bloch state ϕ g −π/d (z) for energy of the seventh miniband bottom, in turn, labelled g (E = 190 meV), are presented. The term d stands for the SL period length. A detailed description of the implementation of the method presented here was published earlier in [20]. The exemplary results for this part of simulation, dedicated to the structure described in [3], are illustrated in Figure 3 where the module (in red) and the phase (in blue) of the selected Bloch state for energy of the seventh miniband bottom, in turn, labelled g (E = 190 meV), are presented.  Figure 4a where the dispersion functions are plotted in the form  Another form of solving Equation (1) is shown in Figure 4a where the dispersion functions are plotted in the form

Another form of solving Equation (1) is shown in
where T ν h is the h th harmonic coefficient for the ν miniband and The term d stands for the SL period length. A detailed description of the implementation of the method presented here was published earlier in [20]. The exemplary results for this part of simulation, dedicated to the structure described in [3], are illustrated in Figure 3 where the module (in red) and the phase (in blue) of the selected Bloch state for energy of the seventh miniband bottom, in turn, labelled g (E = 190 meV), are presented. Another form of solving Equation (1) is shown in Figure 4a where the dispersion functions are plotted in the form  These results illustrate the three lowest minibands calculated for the structure described in [17]. Part b of the Figure 4b shows interesting results, illustrating the conductivity band discontinuity (∆E C ) influence on the miniband position in the energy domain.

The Base of Quantum States Construction
With SLSim, it is possible to construct the base of quantum states, which represent some of the transport properties for specific parameters of the studied structure. For the thermodynamic equilibrium conditions, the base consists of the MLWF [17] While, with the bias voltage applied to SL structure, the program determines the base of the quantum states that consist of the WSF, which are the eigenstates of the Hamiltonian including the electrostatic potential. The methods of the MLWF calculations were described in detail elsewhere [22], while the WSF is determined on the basis of [18]. The results in this area, obtained with SLSim, are shown in Figure 5. These results illustrate the three lowest minibands calculated for the structure described in [17]. Part b of the Figure 4b shows interesting results, illustrating the conductivity band discontinuity (ΔE C ) influence on the miniband position in the energy domain.

The Base of Quantum States Construction
With SLSim, it is possible to construct the base of quantum states, which represent some of the transport properties for specific parameters of the studied structure. For the thermodynamic equilibrium conditions, the base consists of the MLWF [17] While, with the bias voltage applied to SL structure, the program determines the base of the quantum states that consist of the WSF, which are the eigenstates of the Hamiltonian including the electrostatic potential. The methods of the MLWF calculations were described in detail elsewhere [22], while the WSF is determined on the basis of [18]. The results in this area, obtained with SLSim, are shown in Figure 5. The base of MLWF was calculated for the structure described in [3] at the temperature T = 25 K. Based on the above data, it can be concluded that for thermodynamic equilibrium conditions, significant transport in the studied SL does not occur. This is different with the voltage applied, as shown in Figure 6. It shows the WSF base calculated for the same structure for a voltage of 125 mV/period. The distance (in energy) between state 1 (in magenta) and state 2 (in blue) allows the photon emission to occur, while the energy distance between the states 3 (in red) and 4 (in green) opens the possibility of LO-phonon emission. Electron tunnelling through the area of several energy barriers to the next SL period (injection) should be also added to both of the above phenomena, due to the resonance occurring between states 5 (in brown) and 1' (in magenta). The base of MLWF was calculated for the structure described in [3] at the temperature T = 25 K. Based on the above data, it can be concluded that for thermodynamic equilibrium conditions, significant transport in the studied SL does not occur. This is different with the voltage applied, as shown in Figure 6. It shows the WSF base calculated for the same structure for a voltage of 125 mV/period. The distance (in energy) between state 1 (in magenta) and state 2 (in blue) allows the photon emission to occur, while the energy distance between the states 3 (in red) and 4 (in green) opens the possibility of LO-phonon emission. Electron tunnelling through the area of several energy barriers to the next SL period (injection) should be also added to both of the above phenomena, due to the resonance occurring between states 5 (in brown) and 1' (in magenta).
It is arranged in a sequence of so-called cascade electron transport occurring in SL structures and widely used to implement quantum cascade lasers (QCL). However, to accurately determine the scale and type of the quantum phenomena occurring in a working laser, the approach presented here is not sufficient, and it is necessary to apply NEGF formalism, as presented in the next chapter. It is arranged in a sequence of so-called cascade electron transport occurring in SL structures and widely used to implement quantum cascade lasers (QCL). However, to accurately determine the scale and type of the quantum phenomena occurring in a working laser, the approach presented here is not sufficient, and it is necessary to apply NEGF formalism, as presented in the next chapter.

The NEGF Formalism
The majority of information on transport parameters can be found by applying NEGF formalism to model process of the structure. Here, the Dyson and Keldysh equations [12] are solved to determine the retarded G R (E) and correlation G < (E) Green's functions. The Dyson equation can be written as where H is the Hamiltonian of the device, R G is the matrix of the retarded Green's functions and Σ R is the matrix of self-energies. Program initially assumed the self-energies Σ R as the constant diagonal elements iη, where the parameter η is defined with [11] Keldysh equation has the form † R R G G G < < Σ = (8) where G < is the matrix of the Green's correlation functions and Σ < denote the matrix of self-energies.
Initially (in thermodynamic equilibrium) it is the diagonal matrix with elements iη•f n (E), where where the values of E F can be determined in two ways with either the step function, or from an approximation inspired by the Büttiker probes [25,26] applied.
In our approach, a different method for determining G R (E) and G < (E) functions, as well as for the self-energies matrix, is proposed. It is-to the best of our knowledge-different to any of those described to date in superlattice models.

The NEGF Formalism
The majority of information on transport parameters can be found by applying NEGF formalism to model process of the structure. Here, the Dyson and Keldysh equations [12] are solved to determine the retarded G R (E) and correlation G < (E) Green's functions. The Dyson equation can be written as where H is the Hamiltonian of the device, G R is the matrix of the retarded Green's functions and Σ R is the matrix of self-energies. Program initially assumed the self-energies Σ R as the constant diagonal elements iη, where the parameter η is defined with [11] G R (E) = lim Keldysh equation has the form where G < is the matrix of the Green's correlation functions and Σ < denote the matrix of self-energies. Initially (in thermodynamic equilibrium) it is the diagonal matrix with elements iη·f n (E), where where the values of E F can be determined in two ways with either the step function, or from an approximation inspired by the Büttiker probes [25,26] applied. In our approach, a different method for determining G R (E) and G < (E) functions, as well as for the self-energies matrix, is proposed. It is-to the best of our knowledge-different to any of those described to date in superlattice models.
The [26], which describes the exchange of electrons (see Figure 7) between each quantum state of mode b and each electron reservoir of the active device, taking into account the scattered effects considered, was introduced. The definition of the matrix of spectral functions for such an approach takes the form [26] where N indexes all electron reservoirs. It should be mentioned that the diagonal elements of matrix A provide information on the local density of states, which is defined as Electronics 2020, 9, x FOR PEER REVIEW 8 of 17 The quantity ( ) [26], which describes the exchange of electrons (see Figure 7) between each quantum state of mode b and each electron reservoir of the active device, taking into account the scattered effects considered, was introduced. The definition of the matrix of spectral functions for such an approach takes the form [26] where N indexes all electron reservoirs. It should be mentioned that the diagonal elements of matrix A provide information on the local density of states, which is defined as ( ) The transmission coefficient between two electron reservoirs m and n is determined by the relationship [26] [ ] . Trace Knowing the local density of states, and taking into account the scattering effects represented by Γ magnitude, we can determine the local concentration of electrons [26] as where F is a distribution function of electrons, defined by Fermi levels of the Büttiker probes μ N and normalized by the factor k B T. Similarly, the current flowing through a given probe to the quantum state b can be expressed as [26] ( where Fermi levels for the Büttiker probes μ m and μ n are determined from the current continuity  The transmission coefficient between two electron reservoirs m and n is determined by the relationship [26] Knowing the local density of states, and taking into account the scattering effects represented by Γ magnitude, we can determine the local concentration of electrons [26] as where F is a distribution function of electrons, defined by Fermi levels of the Büttiker probes µ N and normalized by the factor k B T. Similarly, the current flowing through a given probe to the quantum state b can be expressed as [26] where Fermi levels for the Büttiker probes µ m and µ n are determined from the current continuity condition and correction in each iteration step of the chemical potentials for each Büttiker probe, according to the relation ∆µ sB = −J −1 I sB (18) until the changes in the µ sB disappear. Having agreed on the values of the chemical potentials of the Büttiker probes, one can determine the self-energies Σ R and Σ < as matrices with diagonal elements in the form iη and iηf Bt (E), respectively, where f Bt . Then, using the Dyson (6) and Keldysh (8) equations, the values of the matrices G R and G < can be calculated. The procedure for reaching these values is ultimately much faster than the one where Born approximation is used, as described in [23]. It is particularly important for simulating the current-voltage characteristics of the device, where each point requires the Green functions' self-consistency. The implementation of above -described method in SLSIM is illustrated in Figure 8.
In Table 2, Fermi levels for the Büttiker probes calculated for the structure described in [24] at the assumed temperature T = 25 K are presented. The impact of voltage on the mixing of quantum states can be observed here, which is one of the most interesting conclusions provided by these calculations.  The calculations with the SLSim are also very effective, as the Hamiltonian matrix in the Equation (6) (19) where H SL and H ξ represent the states of thermodynamic equilibrium and non-equilibrium conditions respectively, while H MF is the part responsible for both scattering of impurities and electron-electron interactions. More details of the method described here were presented in the paper [16], whereas their implementation as part of the SLSim is shown in the Figure 9. As shown, the program allows two options for determining Fermi levels of the quantum states. The first is used to initialise Born's approximation procedure (BAP)) [16] and it is based on the step function. The second one is used in quick procedures for determining GF and it is based on the fixed electron-scattering model (FESM) and Büttiker probes. The use of the BAP allows accurate The calculations with the SLSim are also very effective, as the Hamiltonian matrix in the Equation (6) is small in size, due to the matrix elements being dependent solely on the energy. This Hamiltonian can be represented as a sum where H SL and H ξ represent the states of thermodynamic equilibrium and non-equilibrium conditions respectively, while H MF is the part responsible for both scattering of impurities and electron-electron interactions. More details of the method described here were presented in the paper [16], whereas their implementation as part of the SLSim is shown in the Figure 9. As shown, the program allows two options for determining Fermi levels of the quantum states. The first is used to initialise Born's approximation procedure (BAP)) [16] and it is based on the step function. The second one is used in quick procedures for determining GF and it is based on the fixed electron-scattering model (FESM) and Büttiker probes. The use of the BAP allows accurate calculations of transport parameters, taking into account the given types of electron scattering. This is illustrated by the results shown in Figure 9, where the energy maps of k-resolved occupation functions, calculated on the basis of n(E, k || ) = 2 π ν Im (Tr ([ G < ])), are presented. There are two cases: the first of them illustrates the use of the FESM (in part a) and the second (in part b) takes into account the LO-phonon scattering process calculated with BAP. The simulations are performed for the structure described in [24] at the temperature T = 25 K and a bias voltage equal to 45 mV/period. calculations of transport parameters, taking into account the given types of electron scattering. This is illustrated by the results shown in Figure 9, where the energy maps of k-resolved occupation functions, calculated on the basis of , are presented. There are two cases: the first of them illustrates the use of the FESM (in part a) and the second (in part b) takes into account the LO-phonon scattering process calculated with BAP. The simulations are performed for the structure described in [24] at the temperature T = 25 K and a bias voltage equal to 45 mV/period. The results for the two forms, shown in Figures 6 and 9, combined, are presented in Figure 10 which shows the spatially and energetically resolved electrons' density n(z, E) with the five lowest WSS states.
The simulation was performed for the structure described in [24] with bias voltage equal to 45 mV per period and the assumed temperature T = 25 K. The calculations of the electron density were performed according to [27] ( ) ( ) ( ) ( ), (20) where A is the spatial area in the plane perpendicular to the transport direction, while W α (z) and W β (z) are the states of the Wannier-Stark base.
Unfortunately, a self-consistent solution of the Dyson and the Keldysh equations, as well as self-energy equations in BAP are time-demanding, therefore, when the simulations involve calculating the current-voltage characteristic, the use of the FESM proves to be a more effective approach. Exemplary results obtained with the SLSim are shown in Figure 11. The current-voltage characteristics for the structure described in [24] are presented here. The results for the two forms, shown in Figures 6 and 9, combined, are presented in Figure 10 which shows the spatially and energetically resolved electrons' density n(z, E) with the five lowest WSS states.
The simulation was performed for the structure described in [24] with bias voltage equal to 45 mV per period and the assumed temperature T = 25 K. The calculations of the electron density were performed according to [27] where A is the spatial area in the plane perpendicular to the transport direction, while W α (z) and W β (z) are the states of the Wannier-Stark base. Unfortunately, a self-consistent solution of the Dyson and the Keldysh equations, as well as self-energy equations in BAP are time-demanding, therefore, when the simulations involve calculating the current-voltage characteristic, the use of the FESM proves to be a more effective approach. Exemplary results obtained with the SLSim are shown in Figure 11. The current-voltage characteristics for the structure described in [24] are presented here. Figure 10. The spatially and energetically resolved electrons density n(z, E) with the five lowest WSS states. The simulation was performed with Born's approximation procedure (BAP) for the structure no. 2, described in [24] at the assumed temperature T = 25 K and with the bias voltage equal 45 mV per period.
The total current density in this case can be written as where the matrix [ ] , 0 z H takes into account all direct transitions between states α and β, from the quantum states base currently used, while variable ϑ means the volume of the structure, containing both the number of its periods and the size of the heterojunction surface. The J S current density representing scattering current is calculated as [16] ( where the  The total current density in this case can be written as J = J 0 + J s , where J 0 is the component representing ballistic current, computed according to [18]: where the matrix Ĥ 0 ,ẑ takes into account all direct transitions between states α and β, from the quantum states base currently used, while variable ϑ means the volume of the structure, containing both the number of its periods and the size of the heterojunction surface. The J S current density representing scattering current is calculated as [16] where the dE G < ββ,k (E) αα,k (E) G ad αα,k (E) includes scattering into state α from the state β. The terms in self-energy matrices are set arbitrarily for the scattering type considered currently; in our case they are determined by the value 1 × 10 −4 of the η parameter. It corresponds to the full width at half maximum (∆E fwhm ) of the single peak in the function Im G R ii,k=0 (E) and can be used to measure the decay rate Γ i = ∆E fwhm in scattering process.

Solving the Poisson Equation
With SLSim, it is possible to cover the scattering of impurities and electron-electron interactions alike in the equations by including the H MF part of the device Hamiltonian (see Equation (19)) in the form [16] [ ]  are the creation and annihilation operators, while V MF is the mean field potential, determined by solving the Poisson equation [16] [ ] where ρ d (z) is the structure doping profile, e stands for the electron charge, ε denotes the dielectric constant, while ) (z e ρ is calculated as [16] ( ) ( ) ( ).
Equation (24) is solved with the Newton-Rapson method described in [28]. For this method, Equation (24) (26) where V i is the electrostatic potential in i-th layer of the SL and the dielectric constants ε in the equation are determined in accordance with the relationships Then, the correction is defined in the form Figure 11. The current-voltage characteristics for the structure no. 2. described in []. The total current density J C (plotted in black) is the sum of ballistic current density J 0 (plotted in blue) and scattering current density J 0 (plotted in red). The simulations were performed with the WS base for T = 25 K.

Solving the Poisson Equation
With SLSim, it is possible to cover the scattering of impurities and electron-electron interactions alike in the equations by including the H MF part of the device Hamiltonian (see Equation (19)) in the form [ where m, n and µ, ν indices denote the periods and quantum states of the SL respectively, a ν † m,k,s and a ν n,k,s are the creation and annihilation operators, while V MF is the mean field potential, determined by solving the Poisson equation [16] where ρ d (z) is the structure doping profile, e stands for the electron charge, ε denotes the dielectric constant, while ρ e (z) is calculated as [16] ρ e (z) = 2e Equation (24) is solved with the Newton-Rapson method described in [28]. For this method, Equation (24) is transformed into where V i is the electrostatic potential in i-th layer of the SL and the dielectric constants ε in the equation are determined in accordance with the relationships ε + = (ε i + ε i+1 )/2 and ε − = (ε i + ε i−1 )/2. Then, the correction is defined in the form where j is the index of the grid nodes, while m is the next iteration number. The new potential value after each iteration is therefore Having solved the Poisson equation and calculating the correction (27), the NEGF equations are solved again to take into account the changes in the potential distribution determined by the relation (28), and these actions are repeated until the changes in the potential distribution and charge density disappear (F m i ≈ 0). Sample results for this stage of simulation are illustrated in Figure 12, where the self-consistent charge density n(z) (plotted in blue) and potential E C (plotted in red) simulated for the superlattice structure described in [3] are presented. The calculations are performed for the bias voltage of 215 mV/period at the assumed temperature equal to 25 K. For comparison, the bottom of the conductivity band E C (plotted in black) in the structure prior to starting the charge self-consistent process is also shown in the picture.
where j is the index of the grid nodes, while m is the next iteration number. The new potential value after each iteration is therefore Having solved the Poisson equation and calculating the correction (27), the NEGF equations are solved again to take into account the changes in the potential distribution determined by the relation (28), and these actions are repeated until the changes in the potential distribution and charge density disappear ( Sample results for this stage of simulation are illustrated in Figure 12, where the self-consistent charge density n(z) (plotted in blue) and potential ' C E (plotted in red) simulated for the superlattice structure described in [3] are presented. The calculations are performed for the bias voltage of 215 mV/period at the assumed temperature equal to 25 K. For comparison, the bottom of the conductivity band C E (plotted in black) in the structure prior to starting the charge self-consistent process is also shown in the picture. Figure 12. The self-consistent charge density n(z) (plotted in blue) and potential ' C E (plotted in red) simulated for the structure no. 3. described in [3] with the bias voltage of 215 mV/period and T=25 K.

The Results Presentation
Almost all program blocks illustrated by the algorithm shown in Figure 2 have been programmed in C++ under Microsoft Visual Studio. However, the presentation block for the simulation results was created by using the C # ILNUMERICS library.
The graphic interface allowed us to produce 1D, 2D and 3D chart templates presenting the calculation results, as well as to track the simulations results in real time. The exemplary graphs of Figure 12. The self-consistent charge density n(z) (plotted in blue) and potential E C (plotted in red) simulated for the structure no. 3. described in [3] with the bias voltage of 215 mV/period and T = 25 K.

The Results Presentation
Almost all program blocks illustrated by the algorithm shown in Figure 2 have been programmed in C++ under Microsoft Visual Studio. However, the presentation block for the simulation results was created by using the C # ILNUMERICS library.
The graphic interface allowed us to produce 1D, 2D and 3D chart templates presenting the calculation results, as well as to track the simulations results in real time. The exemplary graphs of the simulation results based on the user-defined templates dedicated for terahertz laser structure [24] are presented in Figure 13. A 1D graph of the summary density of charge on one period structure within energy in the range 0-0.3 eV, and 3D graph of the energy -esolved electron density on the structure period within the energy range 0-0.06 eV, are shown to the left and to the right, respectively. Electronics 2020, 9, x FOR PEER REVIEW 15 of 17 the simulation results based on the user-defined templates dedicated for terahertz laser structure [24] are presented in Figure 13. A 1D graph of the summary density of charge on one period structure within energy in the range 0-0.3 eV, and 3D graph of the energy -esolved electron density on the structure period within the energy range 0-0.06 eV, are shown to the left and to the right, respectively. Figure 13. The exemplary graphs of simulation results based on user-defined templates: a 1D graph (left) and 3D graph (right).

Conclusions
The implementation of the electronic transport simulation methods known from the literature enhanced with our own solutions for Büttiker probes has allowed us to set up a universal and very effective simulator of the heterostructure type I superlattices. When applied, it facilitates a quick assessment of the basic transport parameters of the device under consideration. The program also provides a tool for a thorough analysis of quantum phenomena occurring in superlattices. The simulator has been verified under a number of tests and calculations described in the papers cited therein. The project will be continued for heterostructure type II superlattices, the effects of which will be published.

Conclusions
The implementation of the electronic transport simulation methods known from the literature enhanced with our own solutions for Büttiker probes has allowed us to set up a universal and very effective simulator of the heterostructure type I superlattices. When applied, it facilitates a quick assessment of the basic transport parameters of the device under consideration. The program also provides a tool for a thorough analysis of quantum phenomena occurring in superlattices. The simulator has been verified under a number of tests and calculations described in the papers cited therein. The project will be continued for heterostructure type II superlattices, the effects of which will be published.
Funding: This research received no external funding.