Optimisation of QCL Structures Modelling by Polynomial Approximation

Modelling of quantum cascade laser (QCL) structures, despite a regular progress in the field, still remains a complex task in both analytical and numerical aspects. Computer simulations of such nanodevices require large operating memories and effective algorithms to be applied. Promisingly, by applying semi-analytical polynomial approximation method to computing potential, wave functions and electron charge distribution, accurate results and quick convergence of the self-consistent solution for the Schrödinger and Poisson equations are reachable. Additionally, such an approach makes the respective numerical models competitively effective. For contemporary QCL structures, with quantum wells quite typically forming complex systems, a special approach to determining self energies and coefficients of approximating polynomials is required. Under this paper we have analysed whether the polynomial approximation method can be successfully applied to solving the Schrödinger equation in QCL. A new algorithm for determining self energies has been proposed and a new method has been optimised for the researched structures. The developed solutions have been implemented as a new module for the finite model of the superlattice (FMSL) and tested on the QCL emitting light in the mid-infrared range.


Introduction
QCL structures are nanodevices that have currently grown very popular among scientific teams around the world [1][2][3][4][5]. Their applications likely to be adapted in many important economic branches, such as spectroscopy [6], telecommunications [7], or medicine [8] to name a few, have driven continuous research, with computer simulations playing a major role [9][10][11][12][13][14]. The formalism of non-equilibrium Green's functions has proven to be the most comprehensive approach for it. It can be applied in the form of a real space model (RSM) [15,16], or a model based on the properties of Wannier functions (WFM) [17,18]. The WFM model is faster, and like the RSM, it assumes infinite geometrical sizes of the researched structures. Due to the hardware limitations, simulations with the concerned models are restricted to either one (RSM), or three (WFM) periods of the structure, which, in fact, consists of twenty [19] to nearly two hundred [20] modules. In our previous papers [21][22][23][24] we showed the finite geometrical sizes of the superlattice model to have a significant impact on the simulations results. The finite superlattice model (FMSL) [24] we proposed, due to the semi-analytical approach to the self-consistent solution of the Schrödinger and Poisson equations, not only allowed us to incorporate any number of structural periods, but it was also found as effective as the WFM. The developed FMSL has proven very successful for simulating simple superlattice structures [25] and QCL-THz [20]. Nevertheless, while calculating QCL emitting mid-and far-infrared radiation [26,27], problems aroused. Such devices contain many narrow quantum wells, separated with relatively the source of many photons; hence, such lasers are characterised by high output powers. Figure 1. Illustration of the basic mechanisms for the QCL. The electrons tunnel from the relaxation (injection) region of the previous structure period to the high state (c) of the next superlattice period, followed by a transition to the medium state (b) with a photon emission. In the next step, the electron goes to the low state (a) by emitting a phonon and due to the electric field, it is further transported through the injection area to the high state of the next period (c'), where the sequence of transitions between states is repeated. Namely, we observe the photon transition between states (c') and (b') and then the emission of the phonon after the electron transition from state (b') to (a').
For the correct operation of the device, it is necessary to ensure the population inversion for high quantum states, at a level that allows cascading photon emissions in subsequent periods of the structure. Unfavourable physical phenomena causing the escape of electrons into the energy continuum makes it much more complex. Therefore, the first cascade lasers worked properly only at cryogenic temperatures [29]. The solution to the problem was achieved by appropriately configuring two [30] or three [31] quantum wells, or even complete superlattices [32,33] within the active area of the laser structure. By introducing a very thin quantum well, followed by two wider quantum wells in the active area of the laser, we have managed to significantly increase the efficiency of injecting carriers to the upper level (c) and to reduce the effect of harmful transfer of carriers directly to the lower quantum level (a). All these solutions generate numerous quantum states, often very similar in energies and grouped in very narrow mini-bands. They can be determined by self-consistent solving the Schrödinger and Poisson equations, though it remains quite a difficult task. The problem can be solved in several ways, one of which is described in the following chapter.

The Model of the QCL
The simulations of the QCL structures were carried out in two stages. First, the Schrödinger and Poisson equations were solved self-consistently. Then, the base of quantum states obtained from such calculations was used in while finding the Green functions. On the basis of Green's functions, it was possible to exactly determine the transport parameters characterising the quantum phenomena described in Section 2. The main simulation algorithm is illustrated in Figure 2, where we explicitly show an iterative loop of the self-consistent solution of the Schrödinger-Poisson equation controlled by the value of the min FM δ coefficient, as well as two program blocks responsible for generating the system quantum states base and further determination of Green's func- Figure 1. Illustration of the basic mechanisms for the QCL. The electrons tunnel from the relaxation (injection) region of the previous structure period to the high state (c) of the next superlattice period, followed by a transition to the medium state (b) with a photon emission. In the next step, the electron goes to the low state (a) by emitting a phonon and due to the electric field, it is further transported through the injection area to the high state of the next period (c'), where the sequence of transitions between states is repeated. Namely, we observe the photon transition between states (c') and (b') and then the emission of the phonon after the electron transition from state (b') to (a').
For the correct operation of the device, it is necessary to ensure the population inversion for high quantum states, at a level that allows cascading photon emissions in subsequent periods of the structure. Unfavourable physical phenomena causing the escape of electrons into the energy continuum makes it much more complex. Therefore, the first cascade lasers worked properly only at cryogenic temperatures [29]. The solution to the problem was achieved by appropriately configuring two [30] or three [31] quantum wells, or even complete superlattices [32,33] within the active area of the laser structure. By introducing a very thin quantum well, followed by two wider quantum wells in the active area of the laser, we have managed to significantly increase the efficiency of injecting carriers to the upper level (c) and to reduce the effect of harmful transfer of carriers directly to the lower quantum level (a). All these solutions generate numerous quantum states, often very similar in energies and grouped in very narrow mini-bands. They can be determined by self-consistent solving the Schrödinger and Poisson equations, though it remains quite a difficult task. The problem can be solved in several ways, one of which is described in the following chapter.

The Model of the QCL
The simulations of the QCL structures were carried out in two stages. First, the Schrödinger and Poisson equations were solved self-consistently. Then, the base of quantum states obtained from such calculations was used in while finding the Green functions. On the basis of Green's functions, it was possible to exactly determine the transport parameters characterising the quantum phenomena described in Section 2. The main simulation algorithm is illustrated in Figure 2, where we explicitly show an iterative loop of the selfconsistent solution of the Schrödinger-Poisson equation controlled by the value of the δ FMmin coefficient, as well as two program blocks responsible for generating the system quantum states base and further determination of Green's functions necessary to calculate transport parameters for the structure. The convergence coefficient for the Schrödinger and Poisson equations was defined as where ρ k represents the total charge density of the structure for the kth iteration. The algorithm carries out a loop in which both equations are solved until the charge changes reach the value of δ FMmin .
where ρ k represents the total charge density of the structure for the k th iteration. The algorithm carries out a loop in which both equations are solved until the charge changes reach the value of min FM δ . The numerical model used for the self-consistent solution of both equations is schematically illustrated in Figure 3. We can see there the potential distribution in the direction of electron transport for the QCL structure polarised with the voltage U, along with a set of parameters describing it. In the area of each superlattice layer, a constant value of the effective mass of the electron and a variable potential described by the polynomial function Vj(z) were assumed.  The numerical model used for the self-consistent solution of both equations is schematically illustrated in Figure 3. We can see there the potential distribution in the direction of electron transport for the QCL structure polarised with the voltage U, along with a set of parameters describing it. In the area of each superlattice layer, a constant value of the effective mass of the electron and a variable potential described by the polynomial function V j (z) were assumed.  The numerical model used for the self-consistent solution of both equations is schematically illustrated in Figure 3. We can see there the potential distribution in the direction of electron transport for the QCL structure polarised with the voltage U, along with a set of parameters describing it. In the area of each superlattice layer, a constant value of the effective mass of the electron and a variable potential described by the polynomial function Vj(z) were assumed.   The QCL simulation process begins with solving the Schrödinger equation, as previously reported [24], reduced to a dimensionless form The following dimensionless variables were adopted here as where E, i, and m e j represent the energy and effective mass of the electron in the jth superlattice layer, respectively, while e stands for the elementary charge, and V j (z) for the total potential in the jth layer of the system defined as where V SLj , V Bj , and V S , are the potentials of the superlattice, applied voltage, and unbalanced dopant charge, respectively. Simultaneously, it was assumed that the total potential is represented by a polynomial, which leads to the relationship as below According to the proposed approach, the potential V j (z) is represented by a power series in the form the boundary conditions for Equation (2) were assumed as m e j dΨ j dz z j+1 = m e j+1 dΨ j+1 dz z j+1 (8) while solutions were also sought in the form of a power series where An important problem in the presented approach is to find the number of terms of the series (9), which ensures its convergence to the solutions of the Schrödinger equation (see presented in detail further in Section 5). For the coefficients c j,0 and c j,1 any numerical values can be taken. It allows for obtaining different pairs of independent solutions to Equation (2) in the form Assuming where i = √ −1 and considering that we get After taking into account the boundary conditions we obtain where M j is the transfer matrix for jth layer of the structure in the form for matrix elements described by dependencies and where The transfer matrix M for the whole structure can be written as Assuming zero amplitudes of incident waves from the source and drain sides respectively as we get and finally When solving numerically the Equation (24), the set of self energies of the system is obtained. This provides information on minibands available for electron transport. The task tends to cause problems while simulating specific QCL structures. Therefore, a new computational algorithm (see Figure 4) dedicated to QCL with high energy barriers and narrow quantum wells was developed. The new algorithm is explained in Table 1, where a description for the parameters and procedures is provided. The numerical problems related to the algorithm are described in the following chapters.  Table 1. Basic parameters and numerical procedures of the Schrödinger equation solving algorithm (see Figure 4).

Parameter
Denoted Physical Quantity/Notion Unit 1. Lp The number of the structure periods -2. Vj(z) The total potential in the QCL structure [eV] 3. k Index of the energies currently being considered 4. n Index of the monotonicity vector for the m22(E) function 5. m Index of the self energies 6. Emin The minimum of the assumed energy range The maximum of the assumed energy range Interval of the assumed energy range m2,2(E) function vector 11. NL The number of parts the dEk interval is dived into under MBM procedure -

n KS
The number of terms in series (Equation (9)) under the self energy calculation procedure -

n KD
The number of terms in series (Equation (9)) under the wave function calculation procedure -

n z
The number of the grid points in a single layer of the structure - Transfer matrix procedure   Figure 4).

Lp
The number of the structure periods -

Vj(z)
The total potential in the QCL structure [eV] 3.

n KS
The number of terms in series (Equation (9)) under the self energy calculation procedure - Table 1. Cont.

n KD
The number of terms in series (Equation (9) Optimising the created algorithm provided an important research problem to solve. Still, with parameters defined as no. 11-14 (see Table 1) and following the details described in Sections 4 and 5, we managed to succeed.
As shown in Figure 2, depicting the applied algorithm, the Poisson equation is also solved by the QCL simulation process represented with where ε(z) is the dielectric permittivity, e the elementary charge, and ρ(z) stands for the charge density function, calculated on the basis of the relationships described in our previous paper [24], and by applying the results obtained from solving the Schrödinger equation. Under the Equation (25) both the potential and its derivative at the boundary of each structure layer must be of continuous character, hence The charge density function ρ(z) and the total potential V s (z) in Equation (25) are represented by polynomials in the form [24] The polynomial coefficients were calculated as The last element of the FMSL (see Figure 2) used in the QCL simulation process is the block for determining Green's functions. Its role and operating regime were presented in our previous paper [23]. Encouragingly, unlike the procedures for solving the Schrödinger and Poisson equations, it did not require any additional changes. The obtained results for transport parameters were consistent with the results obtained with alternative models, namely WFM and RSM [22,24] applied. Thus, this fragment of the model is not reported in the paper.

Numerical Analysis of the Function m 2,2 (E)
Two approaches were applied in the process of solving Equation (24). The first one, used also in [24], consisted in discretisation of specific energy range with an appropriate interval and searching for roots of the m 2,2 (E) function by the standard bisection method (BM). The second approach (MBM), presented in this paper, engages the complex function m 2,2 (E) monotony to finding its zeros. Due to large range of values, the tested function was normalised as below It limited the m 2,2 (E)/N function to the range <−1, 1>, which not only significantly facilitated the analysis of the results, but also accelerated numerical calculations.
Application of standard BM to finding the roots of the m 2,2 (E) function, with an appropriately adjusted energy step (dE = 10 −4 ÷ 10 −6 eV), yielded good results for relatively simply superlattice structures, e.g., [24], in terms of obtained results accuracy, and calculations speed. By testing the sign of the function m 2,2 (E), in the selected energy intervals dE k , it was possible to iteratively determine its zeros (see Figure 5a).
It limited the m 2,2 (E)/N function to the range −1, 1, which not only significantly fa cilitated the analysis of the results, but also accelerated numerical calculations.
Application of standard BM to finding the roots of the m 2,2 (E) function, with an ap propriately adjusted energy step (dE = 10 −4 ÷ 10 −6 eV), yielded good results for relatively simply superlattice structures, e.g., [24], in terms of obtained results accuracy, and cal culations speed. By testing the sign of the function m 2,2 (E), in the selected energy inter vals dE k , it was possible to iteratively determine its zeros (see Figure 5a). Applying a similar approach to simulating QCL structures that emit mid and fa infrared radiation failed and resulted in errors, which is illustratively explained in Fig  ure 5b Clearly, at the ends of the energy interval dE k the sign of the function m 2,2 (E) doe not change, so the procedure BM cannot be started. Reducing the energy discretisation to the level as low as dE k = 10 −7 eV not only failed, but also impacted negatively the cal culations efficiency, as n-fold reduction of the dE k means the procedural time for deter mining the roots of the m 2,2 (E) function in the n 2 dimension has been extended.
The problem was overcome by implementing additional numerical procedures which were to discretise the tested energy range with variable intervals, based on the monotonicity analysis of the m 2,2 (E) function. The idea of this approach (MBM) is sche matically illustrated in Figure 6. Applying a similar approach to simulating QCL structures that emit mid and far infrared radiation failed and resulted in errors, which is illustratively explained in Figure 5b Clearly, at the ends of the energy interval dE k the sign of the function m 2,2 (E) does not change, so the procedure BM cannot be started. Reducing the energy discretisation to the level as low as dE k = 10 −7 eV not only failed, but also impacted negatively the calculations efficiency, as n-fold reduction of the dE k means the procedural time for determining the roots of the m 2,2 (E) function in the n 2 dimension has been extended.
The problem was overcome by implementing additional numerical procedures, which were to discretise the tested energy range with variable intervals, based on the monotonicity analysis of the m 2,2 (E) function. The idea of this approach (MBM) is schematically illustrated in Figure 6. culations efficiency, as n-fold reduction of the dE k means the procedural time for deter mining the roots of the m 2,2 (E) function in the n 2 dimension has been extended.
The problem was overcome by implementing additional numerical procedures which were to discretise the tested energy range with variable intervals, based on the monotonicity analysis of the m 2,2 (E) function. The idea of this approach (MBM) is sche matically illustrated in Figure 6. After dividing the interval p1, q1, into NL parts, all the sections are monotonic; hence, th procedure can take the next energy interval p2, q2. The non-monotonicity of the section Lk,2 trig gers recursively its division into NL parts. After dividing the m2,2(E) function into monotonic sec tions, the sign at the ends of each section is checked and the procedure BM is started if signs ar found to differ. After dividing the interval <p 1 , q 1 >, into N L parts, all the sections are monotonic; hence, the procedure can take the next energy interval <p 2 , q 2 >. The non-monotonicity of the section L k,2 triggers recursively its division into N L parts. After dividing the m 2,2 (E) function into monotonic sections, the sign at the ends of each section is checked and the procedure BM is started if signs are found to differ.
Each of the energy intervals dE k is divided into N L sections, and then monotonicity is subjected to tests. For example, in the interval <p 1 , q 1 > the function m 2,2 (E) keeps decreasing for all linear segments, but in the interval <p 2 , q 2 > it is not monotonous within the segment L k,2 . Then, the discretisation procedure is called recursively and the functions on the non-monotonic section are divided into successive N L parts. Once the monotonicity of all newly created sections is tested, the program either invokes recursively the discretisation procedure, or proceeds with the examination of the next energy range dE k . Finally, we obtain a set of energies representing the ends of monotonic sections of the m 2,2 (E) function, which is later used to find its zeros. Then the sign of the m 2,2 (E) function at the ends of its monotonic sections is tested, and if the sign difference occurs, the BM procedure starts (see sections L k,1 i L k,3 ).
The effectiveness of the approach presented here was tested by simulating the structure presented elsewhere [28]. The basic parameters of the simulations and selected results (related to method optimisation) of the MBM are presented in Table 2. It should be added that ∆E denotes the simulated energy range and the parameter T S defines the execution time of the loop for solving the Schrodinger and Poisson equations. Table 2. Basic parameters and selected simulation results for structure presented in [28].

Sym.
No. The results are summarised in Table 3 where the allowed minibands calculated with BM and MBM are provided. Additionally, Figure 7 shows the m 2,2 (E) functions calculated with BM (Sym. 1) and the newly developed MBM. The calculations were performed for 4 and 5 periods of the researched QCL structure ((see Sym. 1) and (Sym. 2), respectively). For all the simulations thermodynamic equilibrium conditions and temperature of T = 100 K were assumed; simulations were performed on a standard PC with Windows 10 and an Intel core i7 processor. As shown in Figure 7, the first detected miniband with the BM applied happens also to be the third consecutive miniband resulting from the MBM approach. It means that the results of Sym. 1 (in blue) have not involved two allowed minibands located in the energy ranges of about 0.06 and 0.1 eV that are vital for electron transport (see states marked as a and b in Figure 1). The picture resulting from with MBM approach looks different, and as shown by the graph of the m 2,2 (E) function obtained in Sym. 2 (in black) it was able to correctly locate all the allowed minibands within the energy range from 0 to ∆E C .

Sym
The research has shown that in our simulations, the N L parameter responsible for dividing the energy range dE k into L k sections (see Figure 6), where the monotonicity of the m 2,2 (E) function is tested, to play a vital role. Its value, if underestimated, may result in failing to detect all the function zeros, as it evidently proved true for Sym. 2, where for N L = 10 in the energy range of 1 mb. (the fragment marked as F in Figure 7a) only 3 out of 5 expected eigenvalues were observed. A very narrow range of this miniband (about 1.3 × 10 −6 eV) and closely located other roots (even within the order of 1 × 10 −8 eV) required increased accuracy of the m 2,2 (E) analysis with regards to its monotonicity. The division of the considered energy range dE k into 20 sections (N L = 20) allowed all the expected self energies to be determined accurately. The research has shown that in our simulations, the N L parameter responsible for dividing the energy range dE k into L k sections (see Figure 6), where the monotonicity of the m 2,2 (E) function is tested, to play a vital role. Its value, if underestimated, may result in failing to detect all the function zeros, as it evidently proved true for Sym. 2, where for N L = 10 in the energy range of 1 mb. (the fragment marked as F in Figure 7a) only 3 out of 5 expected eigenvalues were observed. A very narrow range of this miniband (about 1.3× 10 −6 eV) and closely located other roots (even within the order of 1 × 10 −8 eV) required increased accuracy of the m 2,2 (E) analysis with regards to its monotonicity. The division of the considered energy range dE k into 20 sections (N L = 20) allowed all the expected self energies to be determined accurately. It is worth noting here that applying MBM with fixed dE k for the above calculations is faster and more effective than BM with reduced dE k . For example, the procedure for determining the minibands with N L = 10 lasted 120 s (duration Ts), while with the parameter N L increased to 20, Ts = 1470 s was measured. It is a definite progress when compared to the case where the previous version of the model was used and reducing It is worth noting here that applying MBM with fixed dE k for the above calculations is faster and more effective than BM with reduced dE k . For example, the procedure for determining the minibands with N L = 10 lasted 120 s (duration T s ), while with the parameter N L increased to 20, T s = 1470 s was measured. It is a definite progress when compared to the case where the previous version of the model was used and reducing dE k to 1 × 10 −7 extended the computation time from 80 s (see Table 2) to about 8000 s, without any capability to effectively detect all the allowed minibands in the structure.
As part of the research, the self energy calculations for the structure [3] were also carried out, and their selected results along with the basic simulation parameters are listed in Table 4. The calculations were run under the same supply and environmental conditions as in Sym. 1-Sym. 3. In the subsequent tasks, the parameter N L was increased and the obtained results are shown in Table 4. It turned out that the simulations of the structure [3] needed to increase the accuracy of the calculations related to the structure [28]. It was contributed to the higher potential barriers (∆E c = 0.52 eV) and the related occurrence of two very narrow minibands, instead of one. The first one (1 mb.) lies in the range of 1.3 × 10 −6 eV, while the width of the second one (2 mb.) equals 4.9 × 10 −6 eV. As shown by the results (see Table 4), setting the N L parameter to 20, unlike in Sym. 3, did not provide accurate results. Despite the detection of all 12 allowed minibands, the program failed to calculate a complete set of expected self energies, as for 5 periods of the superlattice, 60 were expected. For N L = 25, the exact values of all self energies were obtained, but as demonstrated, the simulation lasted 5200 s. Table 4. Basic parameters and selected simulations results for structure presented in [3].

Sym.
No. It is also worth mentioning here that calculations where thermodynamic equilibrium conditions are assumed require the highest accuracy, as the smallest self energy difference for the case described above was 9.26 × 10 −9 eV. With voltage applied to the structure, the differences between self energies tend to increase, and consequently, they can be determined much faster.

Convergence of Solutions to the Schrödinger Equation
The proposed solutions to the Schrödinger equation in the form of polynomials raise the question of the number (n k ) of the sequence terms (9) to be applied. A small number of n k means faster calculations, though solutions accuracy suffers. A large number of n k ensures good convergence of the sequence (9) to exact solutions, but it significantly extends computation time. It should be added here that expression (9) is used twice in the process of solving the Schrödinger equation. Primarily it is used to determine the zeros of the m 2,2 (E) function (Equation (24)), and secondly, their application occurs in wave functions calculating procedure for self energies. In the context of optimising computations, it is worth noting that as solving Equation (24) requires greater number of operations, the process lasts much longer than the one for determining the wave functions.
Initial calculations (Sym. 1 to Sym. 3) assumed n k = 100 to ensure very accurate results. It is known, however, that simulations of QCL cover many, often time-consuming, processes (see Figure 2), so it is highly advisable to optimise their duration. It prompted the authors to seek solutions for reducing the sequence length (9) as much as applicable in order to accelerate calculations. The analysis of results obtained for the assumed different values of two parameters gave way to proper approach. Namely, the parameter n kS was responsible for the number of the sequence terms (9) was applied under the process of determining self energies (Equation (24)), while n kD defined the number of such terms for the wave functions calculations.
Applying voltage to the structure also affects the simulation time. For small voltages, self energies within a given miniband differ slightly, particularly for the allowed minibands that are fairly narrow. For such cases detection of quantum states is difficult, as the simulation process is significantly extended. It is also known that for U DS = 0, the base of quantum states is determined, as used in NEGF formalism. Thus, this case seems to be a key factor for the model optimisation process. Taking into consideration all the above, the simulation conditions and two main parameters for the evaluation of their results were successfully defined. The accuracy of the calculated self energies was checked on the basis of the error δ E ν l value defined as where E ν l is the lth value of the calculated self energy within the miniband ν according to sequence (9) of length n kS , while E ν l is the value of the same self energy calculated for n kS = 100. The maximum value Max δ E ν l calculated for l self energies was taken as representative for the miniband ν.
The accuracy of the calculated wave functions was determined in a similar way, i.e., by comparing them to the result obtained at n kD = 100. In this case, the error parameter was defined according to the relationship where Ψ υ l is the lth value of the calculated wave function within the miniband ν according to sequence (9) of length n kD , and Ψ υ l is the value of the same wave function calculated for n kD = 100. The maximum value Max δ Ψ ν l calculated for l self energies was taken as representative for the miniband ν.
Five modules (L p = 5) of QCL structures provided the basis for analysing the obtained results. The assumed number of modules (periods) corresponds to the conclusion reported in our previous paper [23], where it was found to be the minimum size of the FMSL that maximises location of quantum states, and has been found to be an effective base for the NEGF formalism. Another important aspect of the current research was also the tested energy range. We focused solely on the first and eighth minibands (represented by states b and c in Figure 1) due to their fundamental importance for the electron transport in the considered structures. The complexity degree of the respective wave functions that differ significantly has been also considered important for selecting minibands to be analysed. Wave functions representing 8 miniband are much more complicated than those related to 1 miniband. Therefore, most likely the n kD number for miniband c needs to be increased to approximate correctly the wave functions; numerical studies were to answer whether it was necessary.
The semi-analytical approach to QCL modelling allows for calculating the values of the wave functions with high accuracy, without significantly affecting the calculations performance. It is so basically due to the n Z parameter, i.e., the number of calculation points for one layer of the structure.
Selected results of the calculations for the first miniband are illustrated in Figure 8 and listed in Table 5, where the basic parameters of the simulations can also be found. The table contains data from four simulations (Sym. 8-Sym. 11) carried out for different values of parameters related to the length of the sequence (9) (n kS , n kD ). Table 5. Basic parameters and selected results of simulations for structure presented in [28].

Sym.
No. All the plots collectively shown in Figure 8 present the calculated real parts of the wave functions (Re Ψ b ) under thermodynamic equilibrium conditions. The calculations were performed within MBM approach for N L = 20, in the energy range ∆E = 0.06 ÷ 0.07 eV with the basic discretisation step dE = 1 × 10 −6 [eV].

∆E
The results pictured in Figure 8a show that too small number (n kS = n kD = 20) of terms of the sequence (9), despite sufficiently precise discretisation of the energy band (N L = 20), failed to determine all eigenstates in the miniband a. Hence, we observed here only one state with the energy E 1 = 0.06019953 eV instead of the expected five quantum states. The results pictured in Figure 8a show that too small number (n kS = n kD = 20) of terms of the sequence (9), despite sufficiently precise discretisation of the energy band (N L = 20), failed to determine all eigenstates in the miniband a. Hence, we observed here only one state with the energy E 1 = 0.06019953 eV instead of the expected five quantum states. The calculations were performed for the first minibandof the structure consisting of 5 QCL periods reported in [28]. Basic simulation parameters are included in Table 5. Graph  depending on the degree of the approximating polynomial represented by the parameters n kS and n kD . The calculations were performed for the first minibandof the structure consisting of 5 QCL periods reported in [28]. Basic simulation parameters are included in Table 5. Graph (a) shows the wave functions (Re Ψ b ) calculated for the parameters n kS = n kD = 20, graph (b) is the same wave functions calculated for the parameters n kS = n kD = 25, and graph (c) illustrates the above wave functions calculated for the parameters n kS = 25 and n kD = 100. value, but also from relatively wide electrode zone (10 nm was assumed here) in relation to the thickness of other module layers (see Table 2). As the values of the polynomial coefficients are determined at the layers' boundaries, moving away from them when calculating the wave functions, leads to increased inaccuracies.
By increasing the n kS and n kD parameters to 25, we eliminated the failure to detect all the eigenstates, and reduced the error practically to zero (Max δ E b l = 1.5 × 10 −10 %), though good convergence of the wave functions at the ends of the simulated structure remained to be perfected. This is shown by the results of Sym. 10 in Figure 8b, namely, their fragment marked as F, where Max ε Ψ b l = 2.81%. Setting the n kD parameter to the value of 100, as shown in Figure 8c, with the same value of n kS produced correct results for the miniband a and good convergence of the calculated wave functions within the whole structure (Max ε Ψ b l = 0.21%). Correct normalisation of the wave function in the area of the whole structure and the accurate calculation of the related electron charge were possible to be executed in the next stage. It is also worth noting that the parameter selection, proposed in Sym. 11, only slightly increased the computation time from 70 to 78 s, but resulted in nearly six-fold acceleration with respect to the same computations in Sym. 8.
The results for the series (9) convergence with the solutions of the Schrödinger equation in the miniband (c) are shown in Figure 9 and listed with the simulation parameters in Table 6. They include simulations of 5 periods of the QCL [28] under thermodynamic equilibrium with MBM for N L = 20, for the energy range ∆E = 0.23 ÷ 0.24 eV with the basic discretisation step dE = 1 × 10 −6 [eV]. Table 6. Basic parameters and selected results of simulations for structure presented in [28].

Sym.
No. It turned out that setting the parameters n kS and n kD to 20 in Sym. 13, unlike Sym. 9, allowed for detecting all the expected self energies in the miniband under consideration. The maximum error in the self energy value (Max δ E c l ), of only 3.5 × 10 −7 eV occurred here. Regretfully, no good convergence of the function to the assumed zero in the area of the drain, marked as F in the picture, was found and the maximum error (ε Ψ c lmax ) was of 1.55%. This error was minimised to the value of 2.1 × 10 −4 % by setting the parameters n kS and n kD to 25, while the error value Max δ E c l was reduced to 2.4 × 10 −10 eV (see Sym. 14). It all contributed to accurate computation of wave functions and eigenvalues in a very short time of 8 s. It should be added here that the results of Sim. 14 are so close to the results of Sim. 12 that they were not presented in the Figure 9.  (9) with the solutions to the Schrödinger equation (Re Ψ c ) depending on the degree of the approximating polynomial nk. The calculations were performed for the eighth miniband of the structure which consists of 5 QCL periods as described in [28]. Basic simulation parameters are listed in Table 6. Graph (a) shows the wave functions (Re Ψ c ) calculated for the parameters nkS = nkD = 20, graph (b) is the same wave functions calculated for the parameters nkS = nkD = 25.

∆E
It is then reasonable to conclude that despite more complex waveforms of the wave functions in the 8 miniband in relation to 1 miniband, it proved unnecessary to increase the value of the n kD parameter in order determine quickly and accurately all the eigenstates for the tested structure.
The width of the allowed miniband seems to remain the decisive factor with regard to optimising the concerned method. For the concerned case it was relatively large and amounted to 1.89 meV. The subsequent simulations results that were carried out for the structure reported in [3] confirmed such findings. In Table 7 all the basic parameters and selected results of the above calculations carried out within the range of 2 miniband under thermodynamic equilibrium conditions are listed. Table 7. Basic parameters and selected results of simulations for structure presented in [3].

Sym.
No. ceeded the value of the allowed miniband width (4.9 × 10 −6 eV). Hence, the software found only one self energy and the wave function for this energy could not be deter- Figure 9. Convergence of the series (9) with the solutions to the Schrödinger equation (Re Ψ c ) depending on the degree of the approximating polynomial n k . The calculations were performed for the eighth miniband of the structure which consists of 5 QCL periods as described in [28]. Basic simulation parameters are listed in Table 6. Graph (a) shows the wave functions (Re Ψ c ) calculated for the parameters n kS = n kD = 20, graph (b) is the same wave functions calculated for the parameters n kS = n kD = 25.

ΔE
It is then reasonable to conclude that despite more complex waveforms of the wave functions in the 8 miniband in relation to 1 miniband, it proved unnecessary to increase the value of the n kD parameter in order determine quickly and accurately all the eigenstates for the tested structure.
The width of the allowed miniband seems to remain the decisive factor with regard to optimising the concerned method. For the concerned case it was relatively large and amounted to 1.89 meV. The subsequent simulations results that were carried out for the structure reported in [3] confirmed such findings. In Table 7 all the basic parameters and selected results of the above calculations carried out within the range of 2 miniband under thermodynamic equilibrium conditions are listed. Table 7. Basic parameters and selected results of simulations for structure presented in [3].

Sym.
No. The presented results show that for Sym. 16 and Sym. 17 the errors Max δ E b l exceeded the value of the allowed miniband width (4.9 × 10 −6 eV). Hence, the software found only one self energy and the wave function for this energy could not be determined as the sequence (9) was divergent (n.c.). Setting the n kS and n kD parameters to 30 (see Sym. 18) allowed all the expected self energies to be detected and the accuracy of the calculations was increased to the level Max δ E ν l = 1.0 × 10 −7 . A very small error Max ε Ψ b l = 9.3 × 10 −3 required no additional increase for the n kD parameter. In relation to Sym. 15 nearly four times shorter simulation time was achieved.

Conclusions
Contemporary QCL structures emitting mid-and far-infrared radiation typically contain very narrow quantum well systems that tend to generate discrete self energies which are difficult to detect. Hence, numerical models dedicated to such devices need to smartly combine computations that are equally high accurate and efficient. The proposed approach allows a variable discretisation step to be applied with regards to the structure geometrical dimensions and the considered energy range alike. Due to semi-analytical nature of the method, with an appropriate set of parameters for the developed algorithm, simulation process is found optimal, which classifies the FMSL model among highly desirable and effective tools, supportive in designing QCL.