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.

1. 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 high potential barriers. As such, they generate quantum states vital for electronic transport in very narrow minibands, which are hard to detect. Therefore, it was necessary to significantly increase the model accuracy with respect to determining self energies by implementing additional program modules with new numerical algorithms incorporated. The approximation of wave functions resulting from the applied approach was another aspect, as in QCL structures under consideration it can take on very complicated forms. This requires checking and matching the degree of the approximating polynomial while simulating. As the newly developed software significantly extended the simulation time, efforts were focused on optimising the new version of the model.

Under this paper, numerical procedures expanding the FSML application to simulating the transport parameters for QCL emitting mid- and far-infrared radiation have been described and analysed. The solutions to the encountered problems are presented, as well as a significant progress in relation to the previously developed FMSL is shown, in particular, the increased accuracy and efficiency. In order to illustrate the effectiveness of the new model, the simulations results for representative QCL are presented.

As the paper deals with several quite separate issues, we divided it into following sections:

Chapter 2, where operation of QCL with a typical structure emitting radiation in the mid-infrared range is briefly described;

Chapter 3, where we present the numerical QCL model into which a new algorithm for precise self energy search has been incorporated, and its basic procedures are described;

Chapter 4, where validity of implementing such changes into the current QCL model is discussed and supported with relevant calculation results.

And last but not least

Chapter 5, where the convergence of solutions to the Schrödinger equation is analysed, which is of key importance for model optimisation, as confirmed with relevant examples.

The performed studies have finally led to creating an efficient and accurate QCL model dedicated to systems with very narrow quantum wells and high potential barriers for electronic transport.

2. Quantum Cascade Laser

The structure of a quantum cascade laser can be represented as a periodic potential derived from nanometer layers of two different semiconductors arranged alternatively. The application of an electric field in such a system determines electron transport enriched with quantum phenomena, which can be used to emit radiation of different energy. The idea of such a mechanism is schematically illustrated in Figure 1 on the basis of the structure described in [28].

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’).

Each module of the laser structure consists of two parts, namely, the injector, and the active area. Three basic quantum states are available in the active area, namely, a high (c), medium (b), and low (a) state. The electron injected from the previous period of the structure into the state (c) is transferred into the state (b) by emitting a photon of the energy equal to the energy difference between subbands (c) and (b). Then, the carrier reaches the low state (a) and its energy is transferred to the crystal structure as a phonon. The injection area of the next period is designed to carry the electron from the low state (a) of the active region and enable its further transport to the high state (c’) of the next structural module. The presented sequence of transitions means that one electron is the source of many photons; hence, such lasers are characterised by high output powers.

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.

3. 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 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 ρ 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 .

Figure 2

The main algorithm for self-consistent solving the Schrödinger and Poisson equations in the process of QCL simulations.

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 , 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(z) were assumed.

Figure 3

Concepts of polynomial approximation of the potential in the QCL structure. In each superlattice layer, a constant value of the effective mass of the electron and a variable potential described by the polynomial function V(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 represent the energy and effective mass of the electron in the jth superlattice layer, respectively, while e stands for the elementary charge, and V (z) for the total potential in the jth layer of the system defined as where V, V, and V, 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 (z) is represented by a power series in the form the boundary conditions for Equation (2) were assumed as 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,0 and c,1 any numerical values can be taken. It allows for obtaining different pairs of independent solutions to Equation (2) in the form

Assuming where and considering that we get

After taking into account the boundary conditions we obtain where M 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.

Figure 4

The algorithm of the Schrödinger equation solving process under the new approach to determining self energy.

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	[eV]
7. Emax	The maximum of the assumed energy range	[eV]
8. dE	Interval of the assumed energy range	[eV]
9. SE [m]	Self energies vector	[eV]
10. m22 [NL]	m2,2(E) function vector
11. NL	The number of parts the dEk interval is dived into under MBM procedure	-
12. nKS	The number of terms in series (Equation (9)) under the self energy calculation procedure	-
13. nKD	The number of terms in series (Equation (9)) under the wave function calculation procedure	-
14. nz	The number of the grid points in a single layer of the structure	-
Procedure	Denoted Procedure
TM (Ek, nKS)	Transfer matrix procedure
Mon (m22 [ ], NL)	Monotonicity test procedure
Sort (E, M22 [n*NL])	Sorting procedure for the m2,2(E) arguments (ascending order)
Sign (M22[Ek,], M22[Ek+1])	Procedure for checking the sign of a m2,2(E)
BM (M22[Ek,], M22[Ek+1])	Bisection method procedure
WFcalc (SE[m], Lp, nz, nKD)	Wave functions calculation procedure
WFapr ()	Wave functions approximation procedure

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 Section 4 and Section 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(z) in Equation (25) are represented by polynomials in the form [24] and

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.

4. Numerical Analysis of the Function m2,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 m2,2(E) function by the standard bisection method (BM). The second approach (MBM), presented in this paper, engages the complex function m2,2(E) monotony to finding its zeros. Due to large range of values, the tested function was normalised as below

It limited the m2,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 m2,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 m2,2(E), in the selected energy intervals dE, it was possible to iteratively determine its zeros (see Figure 5a).

Figure 5

Schematic illustration of the standard BM engaged in the process of finding the zeros of the m2,2(E) function: (a) an appropriately selected energy discretisation step dEk ensures the sign of the function at the ends of the dEk = segment to be changed and the specific segment to include the root, which can be found in the next iterative step (i + 2); (b) for dEk values that are too large the sign of the m2,2(E) function does not change and the zeros are not found.

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 the sign of the function m2,2(E) does not change, so the procedure BM cannot be started. Reducing the energy discretisation to the level as low as dE = 10−7 eV not only failed, but also impacted negatively the calculations efficiency, as n-fold reduction of the dE means the procedural time for determining the roots of the m2,2(E) function in the n2 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 m2,2(E) function. The idea of this approach (MBM) is schematically illustrated in Figure 6.

Figure 6

A schematic illustration of the new approach (MBM) to finding zeros of the m2,2(E) function. After dividing the interval , into N parts, all the sections are monotonic; hence, the procedure can take the next energy interval . The non-monotonicity of the section L,2 triggers recursively its division into N parts. After dividing the m2,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 is divided into N sections, and then monotonicity is subjected to tests. For example, in the interval the function m2,2(E) keeps decreasing for all linear segments, but in the interval it is not monotonous within the segment L,2. Then, the discretisation procedure is called recursively and the functions on the non-monotonic section are divided into successive N 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. Finally, we obtain a set of energies representing the ends of monotonic sections of the m2,2(E) function, which is later used to find its zeros. Then the sign of the m2,2(E) function at the ends of its monotonic sections is tested, and if the sign difference occurs, the BM procedure starts (see sections Lk,1 i Lk,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 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].

Simulation Parameters								Selected Results
Sym. No.	Sym. Meth.	Lp	Structure Layers [nm]	ΔEc	ΔE	dE	NL(MBM)	Expect. Self En. No.	Calc. Self En. No.	TS[s]
				[meV]
Sym. 1	BM	4	4.6, 1.9, 1.1, 5.4,1.1, 4.8, 2.8, 3.4,1.7, 3.0, 1.8, 2.8,2.0, 3.0, 2.6, 3.0	390	1 ÷ 400	1 × 10−3	-	36	28	80
Sym. 2	MBM	5					10	45	43	120
Sym. 3							20		45	1470

The results are summarised in Table 3 where the allowed minibands calculated with BM and MBM are provided. Additionally, Figure 7 shows the m2,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.

Table 3

Allowed minibands calculated for structure reported in [28] by using the BM method (Sym. 1) and MBM method (Sym. 2 and Sym. 3).

Sym.No.	1 mb.	2 mb.	3 mb.	4 mb.	5 mb.	6 mb.	7 mb.	8 mb.	9 mb.
				[meV]
Sym. 1	125.393 ÷ 125.431	143.827 ÷ 143.961	165.172 ÷ 165.479	187.473 ÷ 187.856	228.578 ÷ 230.212	232.791 ÷ 234.631	340.359 ÷ 341.420
Sym. 2Sym. 3	66.760 ÷ 66.762	108.237 ÷ 108.246	125.390 ÷ 125.431	143.827 ÷ 143.961	165.172 ÷ 165.479	187.473 ÷ 187.856	228.578 ÷ 230212	232.791 ÷ 234.631	340.359 ÷ 341.420

Figure 7

Graphs of the m2,2(E) function obtained in the simulation of the structure reported in [28]: (a) for 4 periods (Sym. 1) and 5 periods (Sym. 2–3) of the superlattice, with BM (Sym. 1) and MBM (Sym. 2–3) applied, with the interval dE = 1−6 eV; (b) calculation results of the m2,2(E) function for the part marked in (a) as F with parameters: N = 10 (Sym. 2) oraz N = 20 (Sym. 3).

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 m2,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.

The research has shown that in our simulations, the N parameter responsible for dividing the energy range dE into L sections (see Figure 6), where the monotonicity of the m2,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 = 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 m2,2(E) analysis with regards to its monotonicity. The division of the considered energy range dE into 20 sections (N = 20) allowed all the expected self energies to be determined accurately.

It is worth noting here that applying MBM with fixed dE for the above calculations is faster and more effective than BM with reduced dE. For example, the procedure for determining the minibands with N = 10 lasted 120 s (duration T), while with the parameter N increased to 20, T = 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 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 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 (ΔEc = 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 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 = 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].

Simulation Parameters								Selected Results
Sym. No.	Sym. Meth.	Lp	Structure Layers[nm]	ΔEc	ΔE	dE	N L	Expect. Self En.	Calc. Self En.	TS[s]
				[meV]				No.
Sym. 4	MBM	5	4.0, 1.9, 0.7, 5.8, 0.9, 5.7, 0.9, 5.0, 2.2, 3.4, 1.4, 3.3, 1.3, 3.2, 1.5, 3.1, 1.9, 3.0, 2.3, 2.9, 2.5, 2.9,	520	1 ÷ 520	1 × 10−3	10	60	51	235
Sym. 5							15		52	373
Sym. 6							20		52	2400
Sym. 7							25		60	5200

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.

5. 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) of the sequence terms (9) to be applied. A small number of n means faster calculations, though solutions accuracy suffers. A large number of n 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 m2,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 = 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 was responsible for the number of the sequence terms (9) was applied under the process of determining self energies (Equation (24)), while n 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 = 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 value defined as where is the lth value of the calculated self energy within the miniband ν according to sequence (9) of length n, while is the value of the same self energy calculated for n = 100. The maximum value 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 = 100. In this case, the error parameter was defined according to the relationship where is the lth value of the calculated wave function within the miniband ν according to sequence (9) of length n, and is the value of the same wave function calculated for n = 100. The maximum value calculated for l self energies was taken as representative for the miniband ν.

Five modules (L = 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 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 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, n).

Figure 8

Illustration of the series (9) convergence with the solutions of the Schrödinger equation (Re Ψ) depending on the degree of the approximating polynomial represented by the parameters n and n. 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 Ψ) calculated for the parameters n = n = 20, graph (b) is the same wave functions calculated for the parameters n = n = 25, and graph (c) illustrates the above wave functions calculated for the parameters n = 25 and n = 100.

Table 5

Basic parameters and selected results of simulations for structure presented in [28].

Simulation Parameters								Selected Results
Sym.No.	ΔE	dE	N L	Lp	n Z	n kS	n kD	Exp. Self En.	Calc. Self En.	Max |δ_Elc|	Max |δ_Ψlc|	TS[s]
	[meV]							No.		[%]
Sym. 8	60 ÷ 70	1 × 10−3	20	5	20	100	100	5	5	-		420
Sym. 9						20	20		1	9.81	72.76	50
Sym. 10						25	25		5	1.5 × 10−7	16.14	70
Sym. 11						25	100		5		0.21	78

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

The results pictured in Figure 8a show that too small number (n = n = 20) of terms of the sequence (9), despite sufficiently precise discretisation of the energy band (N = 20), failed to determine all eigenstates in the miniband a. Hence, we observed here only one state with the energy E1 = 0.06019953 eV instead of the expected five quantum states.

Moreover, the wave function representing the considered state does not converge to the assumed zero in the source and drain regions (see fragments F1 and F2). The maximum error () of convergence of the wave function calculated according to the Formula (33) is 72.76%, whereas the maximum self energy error () for this case equalled 9.81%. Such a large error in the source and drain areas results not only from the high 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 and n parameters to 25, we eliminated the failure to detect all the eigenstates, and reduced the error practically to zero ( = 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 = 2.81%.

Setting the n parameter to the value of 100, as shown in Figure 8c, with the same value of n produced correct results for the miniband a and good convergence of the calculated wave functions within the whole structure ( = 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 = 20, for the energy range ΔE = 0.23 ÷ 0.24 eV with the basic discretisation step dE = 1 × 10−6 [eV].

Figure 9

Convergence of the series (9) with the solutions to the Schrödinger equation (Re Ψ) depending on the degree of the approximating polynomial n. 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 Ψ) calculated for the parameters n = n = 20, graph (b) is the same wave functions calculated for the parameters n = n = 25.

Table 6

Basic parameters and selected results of simulations for structure presented in [28].

Simulation Parameters								Selected Results
Sym. No.	ΔE [eV]	dE [eV]	N L	Lp	n Z	n kS	n kD	Expect. Self En. No.	Calc. Self En. No.	Max |δ_Elc| [%]	Max |δ_Ψlc| [%]	TS[s]
Sym. 12	0.23 ÷ 0.24	1 × 10−6	20	5	20	100	100	5	5	-		41
Sym. 13						20	20		5	3.5 × 10−7	1.55	6
Sym. 14						25	25		5	2.4 × 10−10	2.1 × 10−4	8

It turned out that setting the parameters n and n 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 (), 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 () was of 1.55%. This error was minimised to the value of 2.1 × 10−4% by setting the parameters n and n to 25, while the error value 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.

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 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].

Simulation Parameters									Selected Results
Sym. No.	ΔE [eV]	dE [eV]	N L	Lp		n Z	n kS	n kD	Exp. self en. no.	Calc. self en. no.	Max |δ_Elb| [%]	Max |δ_Ψlb| [%]	TS[s]
Sym. 15	0.1 ÷ 0.101	1 × 10−6	25	5	20		100	100	5	5	-	-	1787
Sym. 16							20	20		1	2.7 × 10−3	n. c.	5
Sym. 17							25	25		1	2.3 × 10−5	n. c.	7
Sym. 18							30	30		5	1.0 × 10−7	9.3 × 10−3	450

The presented results show that for Sym. 16 and Sym. 17 the errors 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 and n 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 = 1.0 × 10−7. A very small error = 9.3 × 10−3 required no additional increase for the n parameter. In relation to Sym. 15 nearly four times shorter simulation time was achieved.

6. 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.