A Polynomial Approximation to Self Consistent Solution for Schrödinger–Poisson Equations in Superlattice Structures

: The paper deals with a new approach to iterative solving the Schrödinger and Poisson equations in the ﬁrst type of semiconductor superlattice. Assumptions of the transfer matrix method are incorporated into the approach, which allows to take into account the potential varying within each single layer of bias voltage superlattice. The key process of the method is to approximate the charge density and wave functions with polynomials. It allows to obtain semi-analytical solutions for the Schrödinger and Poisson equations, which in turn have signiﬁcant impact on the accuracy and speed of superlattice simulations. The presented procedure is also suiﬁhue for ﬁnding eigenstates extended over relatively large superlattice area, and it can be used as an effective pro-gram module for a superlattice ﬁnite model. The obtained quantum states are very similar to the Wannier-Stark functions, and they can serve as the base under non-equilibrium Green’s function formalism (NEGF). Exemplary results for Schrödinger and Poisson solutions for superlattices based on the GaAs/AlGaAs heterostructure are presented to prove all the above.


Introduction
Semiconductor superlattices (SL) are nanostructures whose properties were discovered in the 1970s [1], still the first real devices were produced twenty years later [2].Currently, they are used for building quantum cascade lasers [3,4] which are used in several important industry sectors [5].Due to the high complexity of the physical processes occurring in the SL structures, their computer simulations require either high computing power, or highly effective computational methods.Therefore, over the recent years, many numerical models have been developed for SL structures.In the simplest one [5][6][7][8] a constant potential and the effective mass of electrons in each layer of the structure are assumed.A popular approach to solving the Schrödinger equation in cascade struktures is provided by the rate equation method (RE) [9][10][11][12] where the boundary conditions in the form of hard walls at the edges of the structure composed of three consecutive periods of SL are assumed, and the solutions obtained for the central period are considered to be representative for the whole system.The Monte Carlo method (ME) [13][14][15], based on solving the Boltzman equation, is another approach.Non-equilibrium Green's function formalism (NEGF) [16,17] however, is the approach that offers the most versatile opportunities for simulating superlattices.This approach relies on self-consistent solving of Dyson, Keldysh and selfenergy equations either in the real space (real space method-RSM) [18][19][20], or by means of the Wannier functions (Wannier function method-WFM) [21,22].Under WFM infinite geometrical dimensions of the superlattice, and periodic delocalised Bloch functions over the whole structure are considered.To provide an alternative to such approaches we have built a finite model of the superlattice (FMSL) [23], where the base of the wave functions that is the solution to the Schrödinger equation for the SL model of finite geometrical dimensions is applied.Despite many studies on the semiconductor superlattices properties, conducted both from the empirical and theoretical perspective, many factors that influence the transport parameters of the devices based on such structures still remain unexplored.We believe them to be the ground for the observed discrepancies between the simulation results and measurements performed on real devices.Thus, the pressure to find new tools for studying SL properties under different conditions and fields is still considerable.The influence of finite geometrical sizes on the quantum phenomena occurring in these structures is one of many issues to be studied.Our FMSL [24] where the base of wave functions being the solution to the Schrödinger equation for the superlattice model with a finite geometry is applied, provides quite a suitable tool for such studies.
This paper presents new elements of the FMSL that allows quick studies on the influence of an unbalanced charge, and the voltage applied, on the tested device transport parameters, to be performed.Polynomials used to approximate the charge and potential distribution functions in a SL allow to obtain semi-analytical solutions for the Schrödinger and Poisson equations, where the variable potential in a single structure layer is taken into account without the necessity of dividing it into many sub-areas once the voltage to the system is applied.The presented method is effective in finding eigenstates extended over relatively large superlattice areas.The obtained quantum states are very similar to the Wannier-Stark functions and they can form the basis for the NEGF formalism.In order to confirm such a thesis, exemplary results for Schrödinger and Poisson solutions for superlattices based on GaAs/AlGaAs heterostructure are presented.

Formulation of the Problem
In most cases SL are simply stratified unipolar, n-type nanodevices [25].Figure 1.
shows an example of the AlGaAs-GaAs device presented in [26].Hence, the electrons are subject to the quantum laws, i.e., under the Schrödinger equation with 1D Hamiltonian in the growth d(z) direction within the plane kinetic energy term 2 k 2 /2m e (z), where k = |k| is the magnitude of the in-plane momentum k.In total, it can be written as − 2  2 where Ψ(z) are envelope functions, m e (z) is the effective mass constant within each semiconductor layer, whereas V(z) is a spatially dependent total potential energy which can be written as where V SL (z) represents the superlattice potential, V B (z) is the potential from the applied voltage, and V S (z) denotes the potential derived from impurities and unbalanced charge carriers.
The numerical model of SL presented so far [23,24] have some problems.It consists in the necessity to divide the superlattice layers into many sub-areas with a constant potential (see Figure 2A).In the case of an applied voltage this significantly increases the Hamiltonian device matrix in size.The previous model also does not allow for the consideration of an unbalanced charge.In the present work we propose a solution to these problems by using a polynomial approximation of the potential (see Figure 2B) and then to the charge distribution in the whole SL structure.Many previous calculations (e.g., [27,28]) have shown, that in SL the conduction band is divided into minibands.With the wave functions in the above-mentioned minibands known, it is possible to determine the electron charge distribution within the SL as where n i describes the electrons concentration in the miniband i according to The total charge distribution can be written as where N D is the donor dopant concentration.The charge ρ(z) generates a potential consistent with the Poisson equation where ε is dielectric function.On the border of each layer of the structure, the potential must meet the continuity conditions where z j represents the coordinate of the beginning of the j-th SL layer (see Figure 1), while V is the function of the potential within this layer.The same conditions are fulfilled by the wave function m e j dψ j dz z j+1 = m e j+1 dψ j+1 dz z j+1 (10) In addition, the boundary conditions at the beginning and the end of the structure need to be considered.In our model (see Figure 1) they correspond to the source and the drain, respectively.For the potential function such conditions result from the applied voltage U and it can be written as For the wave function, the conditions correspond to the amplitudes of waves coming from the source and drain, and generally it can be formulated differently [26].All in all, the problem is to determine three coupled fields V, Ψ, and ρ, described by Equations ( 1)-( 6), under boundary conditions ( 7)- (11).

The Main Algorithm of Calculations
The solution of the problem formulated in the previous section requires the fields V, Ψ, ρ, to be self-consistent, which is achievable by calculations with an iterative algorithm applied.The first step is to find the energies and eigenfunctions of the Equation (1) for the V SL and the V B potentials.For V B = 0, the solution of the Equation (1) can be presented as a linear combination of elementary functions, namely plane waves running opposite each other.The standard TMF method is convenient for determining their amplitudes.Incorporating in the next iteration steps the potential of V B , as well as the potential of V S , the solutions to Equation (1) are no longer possible to be presented as elementary functions.Hence, we propose to use solution to the Equation (1) in the form of power series.Such a representation can be obtained analytically, for the potential expressed as a polynomial, or generally, as a power series.By Equation (6) the charge density function takes an analogous form.Therefore, it is suggested to approximate the charge density function with a polynomial, separately for each of the layers, at each step of the calculation procedure, following the solution of Equation (1) and by applying (3)- (5).
The main algorithm of the method under consideration is shown in Figure 3.One of the basic input parameters is the convergence coefficient defined as where ρ k represents the electron charge density for the kth iteration.The algorithm runs a loop for δ > δ FM max .The next programmed operation involves applying the computed wave functions for constructing the quantum states base, to be used in NEGF procedures.
It is worth noting that the self-consistent method for solving the Poisson and Schrödinger equations with the presented algorithm is less complex than the algorithm presented earlier [29], where repeated running of the NEGF blocks have been required that engaged most of the computing power.
The algorithm proposed under the paper significantly speeds up the simulation process and provides an important element of an effective SL model.A detailed description of each of its elements is provided in the following sections.

Solving the Schrödinger Equation
The main task of the procedure described above is to find, separately for each layer of the structure, a general solution to Equation (1) with a given potential in the polynomial form For this purpose, it is advisable to write Equation (1) in dimensionless form, which is obtained by introducing dimensionless variables and For quantities defined his way, Equation (1) takes the form where Solution of Equation ( 16) can be also represented as a power series Substituting ( 17) and ( 18) to ( 16) a recursive relationship of the coefficients c j,n is obtained The coefficients c j,0 and c j,1 can be assumed arbitrarily, instead of adopting specific numerical values, which allows to obtain different pairs of independent solutions to Equation (16).For further analysis, a pair of functions F I j (u), and F II j (u) defined by dependencies ( 18), ( 19) and where i = √ −1 are applied.By this choice, compliance with the solutions representing counter-rotating waves is achieved for . The general solution of Equation ( 16) is presented as Under continuity conditions ( 9), (10) we get the following dependencies where and The Equations ( 23) and ( 24) can be transformed with the matrix notation where is the state-transition matrix for jlayer of the structure.By using (27) we get where is the state-transition matrix for whole structure.In the Equation ( 29) the constants C I 0 and C II L+1 are the amplitudes of incident waves from the source and drain respectively.According to the assumption C I 0 = 0, C II L+1 = 0, by substituting to (30) we arrive at where m 1,2 and m 2,2 are the elements of the matrix M. By solving Equation ( 32), numerically we obtain the set of eigenenergies of the system under analysis.

Solving the Poisson Equation
According to the assumptions, the right side of the Equation ( 5) is approximated in the j-th superlattice layer by an N-degree polynomial The approximation procedure consists in comparing the charge density function calculated in the previous iteration step with the function described by Equation ( 33) at N + 1 selected points (nodes) of each structure layer.This leads to the algebraic systems of equations of N + 1 degree, by solving which, the values of the coefficients a j,k are obtained.Substituting (33) to (6), and following double integration, we get Comparing (33) with (12) we get Undetermined coefficient values b j,0 , and b j,1 under Equation (34) can be found by applying the boundary conditions ( 7), (8) Equations ( 36) and (37) can be written in matrix form where Equation ( 38) allows coefficients b j,0 and b j,1 in the subsequent layers to be calculated recursively and along with (34) it is possible to determine the full potential distribution within the superlattice structure.The coefficient values for the first layer (b 1,0 and b 1,1 ) should be determined from the condition (11).For this purpose, it is necessary to link them with the coefficients of the last layer (b L,0 , b L,1 ).Using (38) we get where with Matrix P in (42) is written: where Using (11), (34) and (38) we get b 1,0 = 0, (47 As we can see, the proposed method for determining the potential distribution within the superlattice consist mainly in performing multiplying and adding operations on 2-rank matrices, which highly speeds up calculations.

NEGF Modelling of Superlattice Structure
The FMSL under the paper incorporated the NEGF formalism into calculating transport parameters.As a part of the NEGF procedure, the Dyson and Keldysh equations are solved to determine the retarted Green's (G R ) and correlation (G < ) functions, respectively.the model described here has an implemented procedure for calculating Green's function based on the WFM method [29] where we solve Dyson equation in the form 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.The Keldysh equation takes the form where G < is the matrix of non-equilibrium Green functions, also called density or correlation functions.The elements of the self energy matrices (Σ < and Σ R ) are found in the procedure of solving self-consistent solution equations together with the matrix elements G < and G R , respectively for the considered electron scattering mechanisms.This allows to determine transport parameters of the tested device, such as the density of state (DOS) [23].
and electron density Transport parameters under this study were calculated with a module into which the Büttiker probes had been implemented [25].It is a well known approach to solving Equations ( 49) and (50), the essence of which is the exchange of electrons between each quantum state of the active device and the virtual electrode attached to it (a Büttiker probe) containing the electron reservoir.As the result, self-consistent Fermi levels for the quantum and electrochemical state of the electrode are attained.It allows to determine the occupation of states in a wide energy range.It is also worth noting that such a method for determining Green's function is very effective, as shown in our earlier paper [28].
Under another research [29], we found that the calculation speed can be significantly reduced when the approach with the WFM method is applied, according to which the Hamiltonian matrix is written as a sum, H = H 0 + H R with the matrix H 0 representing the sum of states for the non-polarized (H SL ) and polarized (H ξ ) superlattice, and the matrix H R containing all electron scattering effects.Detailed definitions of all of the above Hamiltonian matrices are given in our previous papers [29,30] whereas this article has been focused on one of them only, related to the applied voltage, which appears in the form where e is the electron charge, ξ is the electric field, n stands for the number of the considered superlattice period of length-d, k is the wave vector, and the parameters > a ν † n,k and > a ν n,k are state operators, namely creation and annihilation of the ν state, respectively.The variable R µν l included in the definition (51) calculated as is called overlapping states integral.The interaction between the minibands (indices µ, ν) and the superlattice periods (index l) are taken into account by this relation.Under IMS model [28] the W (z) functions correspond to the Wannier and Wannier-Stark basis for unpolarised and polarised superlattice structures, respectively, whereas under the FMS model they correspond to the constructed quantum states described elsewhere [23].Under this paper, we use the base of eigenstates that was created by solving the Schrödinger equation on the basis of the procedure presented elsewhere [23].

Simulations Results
The correctness and effectiveness of the developed numerical method has been checked by simulating typical superlattice structures whose parameters are widely available elsewhere [26,28] that for the purpose of this paper are called structure 1 and structure 2, respectively.They are presented in Table 1.

Materials of Heterostructures
The An important part of numerical research was concerned with assessing the impact of the SL geometrical sizes on the quantum states properties.The results of the research in this area are illustrated in Figure 4.The energy eigenvalues distributions calculated with the FMSL are clearly dependent on the number of the superlattice periods considered.Picture (A) presents the results for a simple symmetrical superlattice that consists of one well and one barrier in the period (structure 1), while picture (B) displays the eigenvalues of energy for a typical terahertz laser structure (structure 2).The results show a general tendency for the energy eigenvalues to group within the minibands.The minibands boundaries are defined by the number of SL periods.This is very clearly seen in the graphs presented in Figure 4A, thought it is less clearly pronounced in Figure 4B, due to the very narrow minibands for this SL.In Figure 4. it can also be observed that the greater number of the periods considered in the calculations, the more precisely the minibands edges are determined.It indicates the necessity of considering the greatest possible number of SL periods in simulations.Regretfully, it practically brings very high load to calculating machines, as well as decreases simulations efficiency.
Therefore, a question arises on the SL periods number to be taken into account in simulations that would give the exact results for the miniband distribution, without the efficiency loss for calculating SL parameters.A comparison of the presented results with the results obtained by means of IMSL may provide an answer.To compare we used the coefficient defined as where υ δ -is the value of the parameter for the υ -miniband calculated using the It is also worth noting that the eigenstates located near the upper edge of the minibands exceed the energy barrier levels (∆E C ) already for 4 periods for structure 1, and 8 periods for structure 2. When these states are occupied by electrons, it may have a significant impact on transport.
In Figure 4. it can also be observed that the greater number of the periods considered in the calculations, the more precisely the minibands edges are determined.It indicates the necessity of considering the greatest possible number of SL periods in simulations.Regretfully, it practically brings very high load to calculating machines, as well as decreases simulations efficiency.
Therefore, a question arises on the SL periods number to be taken into account in simulations that would give the exact results for the miniband distribution, without the efficiency loss for calculating SL parameters.A comparison of the presented results with the results obtained by means of IMSL may provide an answer.To compare we used the coefficient defined as where δ υ I M -is the value of the parameter for the υ -miniband calculated using the IMSL, while δ υ FM it is the value of the same parameter determined by using of the FMSL.From important parameters that characterize superlattices the minibands width (∆E) and the location of the miniband centers energy (E center ) were selected to compare both models.Then the values were determined for each miniband of the considered structure, and their dependence on the number of SL periods with the FSML model, and compared to the results obtained with IMSL.The results are illustrated in Figure 5.Both models were found to converge nicely for simulations with 40 SL periods taken into account under FMSL.The differences increased to less than 1%, and to 0.1% for the width of the minibands (∆E) and the position of their centers (E center ), respectively.It was observed for both simulated structures.It is worth noting here that the obtained results are also consistent with the results published in [29].In addition practical implementations for structure 1 contain several hundred periods, whereas for structure 2 the number of periods to be considered is equal to 190 [26].Hence, any errors due to assumed infinite sizes in IMSL seem negligible.On the other hand, it can be also argued that slight differences observed in the results are encouraging with regards to applying the FMSL to a greater extent than it has been practiced so far.It is clear that the developed method allows to engage the structure bias voltage directly in the process of solving the Schrödinger and Poisson equations.The calculations showing the results are presented in Figure 6A shows the bias voltage influence on the eigenenergy distribution for the structure 2, calculated by means of the FMSL.The effect of so-called mixing eigenstates from different minibands can be observed.For example, from the voltage of 5 mV/period, the energy for the highest state of miniband b 2 is greater than the lowest state of miniband c 2 (mix bc).From the voltage of 6 mV/period, it can be seen that the highest state of the d 2 miniband occupies a higher energy level than the lowest state of the e 2 miniband (mix de).For another case, the energy of the highest state of miniband a 2 is greater than for the lowest state of miniband b 2 (mix ab) already from the voltage of 7 mV/period, and from the voltage of 9 mV/period, the energy for the highest state of miniband c 2 is greater than the lowest level of miniband d 2 (mix ab).All the cases can result in the interband transitions combined with tunneling effects through the potential barriers.These are the electronic transport features characteristic for this type of SL structures.Potential occurrences of such quantum effects are illustrated with Figure 6B, where we can see the wave functions (modules squared) for the bias voltage of 10 mV/period applied to the structure.
Mixed quantum states with similar eigenvalues, for example d iii 2 and e ii 2 , drawn in green and pink, respectively, located in adjacent quantum wells separated by narrow potential barriers are found there.It enables the electron tunnel transport in the local area of the superlattice, even as far as to the next period, to the energy level from the higher mini-band.This is a process known as charge carrier injection.In the next stage, electrons from a higher energy level due to local inversion of occupations can be subject to a coherent transition to the average level (for example d ii 2 → c ii 2 ) emitting a photon with an energy characteristic for the simulated structure (here 14 meV).Vertical transitions incoherent between quantum states emitting LO phonons with the energy of 37.5 meV, such as transitions b ii 2 → a ii 2 , are also possible, with appropriately biased voltage.Such processes result in a electron cascade flux transporting the charge along the superlattice structure, and they provide the basis for cascade lasers operation.However, the precise conditions for such transport can be established only following the NEGF formalism.
In order to perform such calculations, a base of quantum states adequate to the current bias voltage should be constructed and applied to creating the tested device Hamiltonian matrix (51), which comprises the input data of the program block calculating Green's functions.Calculations results for several various quantum states bases are known [31,32] and they can serve as the bases for positional quantum states, the quantum states Wannier functions, or the Wannier -Stark functions.So far all the models developed by our team have applied the Wannier and Wannier-Stark base for the IMSL [29], and the base of quantum states defined for the FMSL [23].The disadvantage of the latter is that the quantum states representing individual minibands are not localized within one SL period.It implies the necessity of including at least five structure periods in the tested device Hamiltonian matrix, while in the IMSL only three structure periods are required for the same purpose.In the approach under the paper, all the obtained quantum states are very well localized within a single SL period, as is the case with the base of quantum states of the IMSL.Hence, the size of the Hamiltonian matrix can be limited to 3 SL periods, as in the IMSL.
With the quantum states base shown in Figure 6 and the size of the Hamiltonian matrix limited to three SL periods, the transport parameters for structure 2 are simulated and selected examples of the results are shown in Figure 7.   b 2 ), so electronic transitions emitting photons with an energy of 14 meV may occur.Additionally, transitions in the k-vector space emitting optical phonons with an Maps of DOS (N(E, k || )) and occupation by electrons (n(E, k || )) functions in the energy space of the wave vector k are also presented.The calculations run at the assumed temperature T = 77 K.The graphs (A) and (B) were made for the bias voltage of U = 10 mV/period, while the graphs (C) and (D) are related to the voltage U = 50 mV/period.By analyzing the results, it can be concluded that for a small voltage applied to the structure (10 mV/period) we observe a high density of states in the energy range corresponding to their eigen energies (see Figure 6B).Electrons occupy only the lowest states (a 2 and b 2 ) and there is no visible inversion of occupations opening the possibility for photons emission (transitions d 2 → c 2 ), or phonons (transitions d 2 → a 2 ).The situation is different for the superlattice biased with voltage of 50 mV/period (see Figure 7D).The higher energy states (e.g., d ii 2 ) show then a greater occupation than the mid-states (b ii 2 ), so electronic transitions emitting photons with an energy of 14 meV may occur.Additionally, transitions in the k-vector space emitting optical phonons with an energy of 37.5 meV are also likely to occur.More detailed research with NEGF applied in this area was published earlier [23,29].This article has been aimed at testing the correctness of a newly developed computational module, whose main task is to effectively find self-consistent potential and charge in a biased voltage SL structure.For this purpose, simulations for typical SL structures were carried out, the results of which are presented in Figure 8A shows the results of a self-consistent solution for the Schrödinger and Poisson equations with regards to potential and total charge distributions in the structure 2 with a bias voltage at 50 mV/period; the structure has been doped at N D = 4 × 10 24 m −3 .Graph (C) shows convergence of the method for different doping levels, i.e., δ FM coefficient (see Equation ( 12)) versus the iteration number dependencies.The convergence of the method was tested in a relatively wide range of impurity concentrations, namely from N D = 2 × 10 23 m −3 to N D = 6 × 10 24 m −3 .The results prove the method to be quickly convergent and to require a maximum of 30 iterations, with the minimum charge change of δ FMmax < 0.1% assumed.Simulations for the same structure with the input parameter defined with the Newton-Rapson method [33] were performed for comparison purposes.The results are represented by charts (B) and (D).Chart (B) presents the potential energy E C distributions and the total charge ρ(z) for the polarized structure, doped at the level of N D = 1.9 × 10 23 m −3 , while graph (D) displays the dependence between δ IM coefficient and the number of iterations.The convergence coefficient δ IM is calculated from where k is the iteration number, δV max stands for the input parameter assumed at 10 −4 V, and δV denotes the potential correction calculated according to the procedure presented in [33].By analyzing the results shown in Figure 8A in relation to the graphs presented in Figure 8B correctness of the new method can be considered confirmed, as similar results with regards to both total charge density distribution and its influence on the E C potential were obtained.With regards to the results shown in graphs (C) it should be noted that the very high doping N D = 6 × 10 24 m −3 has not prevented us from achieving good convergence.However, it is not the case for the Newton-Rapson method.As seen from chart (D) a relatively low level of doping, namely N D = 1.9 × 10 23 m −3 , requires 34 iterations for the charge and potential within the structure to be elf-consistent under the conditions specified in Figure 8, and for very high doping levels convergence problems may occur.Similar cases have been reported [34] and PID blocks were recommended for the method algorithm in order to obtain stable and self-balanced procedure for the charge and potential, which in turn is associated with a significant extension of the simulation time.
The linear relationships on Figure 8C suggest that the convergence of the developed method is exponential (ordinate axis is on a logarithmic scale).This is its great advantage, because it means fast convergence as well as allows to predict simulation time.

Conclusions
The method presented here introduces a significant change in the approach to solving the Schrödinger and Poisson equations in I type superlattices, which-to the best of our knowledge, has not been used in any model for such structures so far.
The polynomial approximation of the charge density function allows to obtain analitycal solutions of the Poisson and Schrodinger equations, which significantly accelerates the procedure for obtaining self-consistent fields in tested system.
As the method was tested on relatively simple structures, with maximum 8 layers in the period, the study of more complex structures remains to be carried out.It requires, however, a better adjustment of polynomial approximation coefficients to more complex wave functions.This topic is currently under research.

Figure 1 .
Figure 1.A fragment of AlGaAs-GaAs superlatice and illustration of periodic potential (with period d and number of the semiconductor layers L), which it is created in the material as a result of the modulation arrangement of the semiconductor layers.The layer sequence of one period in nanometers is: 7.8, 2.4, 6.4, 3.8, 14.8, 2.4, 9.4, 5.4 nm.AlGaAs layers are denoted in bold and the underlined layers is n doped.

Figure 2 .
Figure 2. Concepts of approximation of the potential in the SL: (A) constant potential in each sub-layer of the SL, (B) polynomial approximation of the potential in whole SL.

Figure 3 .
Figure 3.The main algorithm of self-consistent solving of Schrödinger and Poisson equations in the first type of superlattice structures for FMSL.

Figure 4 .
Figure 4. Energy eigenvalues distribution dependence on the number of SL periods calculated with FMSL for structure 1 and structure 2 are presented in parts (A) and (B), respectively.The height of the energy barriers (ΔE C ) is marked blue.

Figure 4 .
Figure 4. Energy eigenvalues distribution dependence on the number of SL periods calculated with FMSL for structure 1 and structure 2 are presented in parts (A) and (B), respectively.The height of the energy barriers (∆E C ) is marked blue.

Figure 5 .
Figure 5.The convergence for IMSL and FMSL illustrated by exemplary minibands simulation for two typical SL structures.The dashed line corresponds to the results for structure 1, while the dependencies concerned with structure 2 are plotted with solid ones.Part (A,B) show the convergence of the width of minibands (∆E), and the energy positions at the centers of these minibands (E center ), respectively.

Energies 2021 ,
14, x FOR PEER REVIEW 15 of 18 With the quantum states base shown in Figure 6 and the size of the Hamiltonian matrix limited to three SL periods, the transport parameters for structure 2 are simulated and selected examples of the results are shown in .

Figure 7 .
Figure 7. Simulations for the DOS ( ) , ( || k E N ) and the electron density ( ) , ( || k E n ) functions in the energy space of the wave vector k.The calculations were performed for structure 2 with NEGF and the quantum state base of the FMSL, for temperature T = 77 K. Graphs (A) and (B) correspond to voltage U = 10 mV per period, whereas graphs (C) and (D) refer to the voltage U = 50 mV per period.Maps of DOS ( ) , ( || k E N ) and occupation by electrons ( ) , ( || k E n ) functions in the energy space of the wave vector k are also presented.The calculations run at the assumed temperature T = 77 K.The graphs (A) and (B) were made for the bias voltage of U = 10 mV/period, while the graphs (C) and (D) are related to the voltage U = 50 mV/period.By analyzing the results, it can be concluded that for a small voltage applied to the structure (10 mV/period) we observe a high density of states in the energy range corresponding to their eigen energies (see Figure 6B).Electrons occupy only the lowest states ( 2 a and 2 b ) and there is no visible inversion of occupations opening the possibility for photons emission (transitions 2 2 c d → ), or phonons (transitions 2 2 a d → ).The situation is different for the superlattice biased with voltage of 50 mV/period (see D).The higher energy states (e.g.ii d 2 ) show then a greater occupation than the mid-states ( iib2 ), so electronic transitions emitting photons with an energy of 14 meV may occur.Additionally, transitions in the k-vector space emitting optical phonons with an

Figure 7 .
Figure 7. Simulations for the DOS (N(E, k || )) and the electron density (n(E, k || )) functions in the energy space of the wave vector k.The calculations were performed for structure 2 with NEGF and the quantum state base of the FMSL, for temperature T = 77 K. Graphs (A,B) correspond to voltage U = 10 mV per period, whereas graphs (C,D) refer to the voltage U = 50 mV per period.

Figure 8 .
Figure 8. Results of self-consistent solution to Schrödinger and Poisson equations, i.e., distribution of total charge-ρ and potential energy-E C ) for structure 2 with bias voltage at U = 50 mV per period and the doping at the level of N D = 4 × 10 24 m −3 obtained with FMSL (part (A)), as well as doping at the level of N D = 1.9 × 10 23 m −3 by means of IMSL (part (B)).Graph (C) shows the dependence of the convergence coefficient-δ FM on the iteration number.The individual series in the chart refer to a different doping of the SL structure.Graph (D) displays the dependence of the convergence coefficient δ IM on the iteration number for the structure doping at the level of N D = 1.9 × 10 23 m −3 .All the computations were performed for the temperature T = 77 K.