Modeling Energy Bands in Type II Superlattices

: We present a rigorous model for the overall band structure calculation using the perturbative k · p approach for arbitrary layered cubic zincblende semiconductor nanostructures. This approach, ﬁrst pioneered by Kohn and Luttinger, is faster than atomistic ab initio approaches and provides sufﬁciently accurate information for optoelectronic processes near high symmetry points in semiconductor crystals. k · p Hamiltonians are discretized and diagonalized using a ﬁnite element method (FEM) with smoothed mesh near interface edges and different high order Lagrange/Hermite basis functions, hence enabling accurate determination of bound states and related quantities with a small number of elements. Such properties make the model more efﬁcient than other numerical models that are usually used. Moreover, an energy-dependent effective mass non-parabolic model suitable for large gap materials is also included, which offers fast and reasonably accurate results without the need to solve the full multi-band Hamiltonian. Finally, the tools are validated on three semiconductor nanostructures: (1) the bound energies of a ﬁnite quantum well using the energy-dependent effective mass non-parabolic model; (2) the InAs bulk band structure; and (3) the electronic band structure for the absorber region of photodetectors based on a type-II InAs/GaSb superlattice at room temperature. The tools are shown to work on simple and sophisticated designs and the results show very good agreement with recently published experimental works.


Introduction
In recent years, developments in crystal growth techniques and in semiconductor device technology has resulted in the fabrication of complicated nano-heterostructures with interesting physical phenomena and promising device applications. Thus, new methods and tools are required for the simulation and analysis of such nanostructures [1]. Among these nanostructures, Antimony-based type-II superlattices (T2SLs) have gained a lot of attention in the field of infrared detectors and lasers due to their high operating temperature. T2SLs combined with unipolar barrier architectures [2,3] have been recently demonstrated to operate at room temperature in infrared photodetectors [4] and interband cascade lasers [5]. Another interesting property of T2SLs is the ability to tune an effective band gap over a wide range by varying layer thicknesses instead of mole compositions [6,7]. The accurate determination of the overall band structure and wave functions of nanostructures is very complicated to tackle analytically except for a few special cases. Therefore, numerical calculations are required to handle complex designs. Among the difficulties of calculating the full band structure for heterostructures is the sensitivity of modeling tools and techniques to the input parameters, namely effective masses, band gap energies, and valence band offset (VBO), especially for structures made of narrow gap materials such as InAs and InSb [7]. Importantly, the exact VBO values are crucial for modeling heterostructures, yet they are very sensitive to the heterointerface imperfections and even a small shift leads to erroneous results. The choice of the modeling technique often comprises a trade-off between accuracy and computational load. Generally, the approaches used for the overall band structure calculation fall into two categories: first principle calculations (i.e., ab initio) [8,9] and empirical methods [10][11][12]. Ab initio calculations such as density functional theory and the quantum Monte Carlo start from atomistic values where information about atomic orbitals and atoms' positions are required without the need of any fitting and/or empirical parameters [8,9]. Although these methods yield highly accurate results, they require considerably long execution time and hence they are limited to systems of a small number of atoms. In contrast, empirical methods are usually used for band structures calculation of large semiconductor heterostructures where some properties are approximated using experimentally fitted parameters. These fitting parameters give a tremendous reduction to the computational cost. One of the most used empirical methods for the band structure calculation is the k · p method, which treats heterostructures as a stack of bulk materials. The method is convenient to model large complex designs due to its highly computational efficiency and good accuracy, despite the hidden atomic scene. The aim of this work is to develop rigorous modeling tools for calculating the band structure of nanostructures using the k · p approach. Since this method is sensitive to variations of input parameters across the heterointerface, the finite element method (FEM) is chosen as a numerical technique to discretize the k · p Hamiltonians. The finite element method has been applied to many quantum mechanical problems ranging from simple isolated systems such as the hydrogen atom [13][14][15] and the single quantum well [16,17] to more complex problems such as the periodic potential [18] and shown to give very precise numerical results in calculating the energies of the system. However, due to the complexity of the FEM and the complicated formulation of the variational functional for a given multiband k · p Hamiltonian, finite difference method (FDM) is usually favored over FEM. Despite FEM accuracy, its application or use remains limited. In FEM, a non-uniform mesh as well as high order Lagrange/Hermit basis functions can be utilized [19], hence the band structure and wavefunctions for heterostructures with layers of arbitrary thicknesses can be calculated accurately with a smoothed mesh the near heterointerface edges. In this work, we first choose the Luttinger-Kohn Hamiltonian as a case study and present its formulation using k · p theory taking into account strain effects and the coupling between conduction and valence bands. This formulation is adapted to heterostructures by using the envelope function approximation (EFA). Then, the finite element discretization scheme for the k · p Hamiltonian is presented by applying the variational principle [20] to the multiband Hamiltonian, which leads to an equivalent eigenvalue matrix representation of the problem. Next, the accuracy of the tools using energy-dependent effective mass non-parabolic model is investigated and compared to published results on a structure of a single quantum wells using FEM, FDM, and the Transfer Matrix Method (TMM). Then, the tools are applied to calculate the bulk dispersion of InAs material under compressive strain. The tools are also applied to calculate miniband dispersions and absorption coefficient of the active region of an infrared photodetector based on type-II InAs/GaSb superlattices operating at room temperature.

k · p Band Structure Theory
Based on the theory of Luttinger-Kohn [21] and Bir-Pikus [22], the 8 × 8 k · p Hamiltonian for bulk zincblende crystal, which describes the interaction between the conduction and the valence bands (heavy-hole (HH), the light-hole (LH), and the spin-orbit split-off (SO)) including strain effects reads: where The shape of the Hamiltonian usually depends on the states and on the order of the basis states used for the construction. The above Hamiltonian is formulated using the following basis states where it is assumed that the growth z-axis of the crystal coincides with the crystallographic axis [001]: These states are formed by an appropriate linear combinations of atomic p-like zone center basis states |X ↑ , |Y ↑ , |Z ↑ and their spin counterpart |X ↓ , |Y ↓ , |Z ↓ , respectively, as described by Kane [23].
In addition, γ 1 , γ 2 , and γ 3 are Kohn-Luttinger parameters; E c and E v stand for the unstrained conduction and valence band edges, respectively; k x , k y , and k z are the wave vector components; k 2 t = k 2 x + k 2 y ; m * e is the conduction-electron effective mass and m * 0 is the free electron mass; P cv is the Kane's interband momentum-matrix element parameter, which is also related to the Kane's energy E p ; ∆ is the spin-orbit splitting energy; a c , a v , and b are the Bir-Pikus hydrostatic and shear deformation potentials, respectively; and xx , yy , and zz are the nonzero components of the strain tensor in the case of a biaxial Strain. This case covers most of the important strained systems, namely a bulk semiconductor under an external uniaxial stress and a semiconductor strained-layer pseudomorphically grown on a (001) substrate. These components of the strain tensor arise from two strain contributions. Firstly, intrinsic due to lattice constant mismatch defined by: where a sub and a are the lattice constant of the substrate and of any layer material within the heterostructures, respectively. C 11 and C 12 are the elastic stiffness constants. Secondly, extrinsic due to external applied stress, which, for the case of an external uniaxial stress [100] : xx = S 11 T, yy = zz = S 12 T. [110] [001] : xx = yy S 12 T, zz = S 11 T.
where S 11 , S 12 , and S 44 are compliance constants related to the elastic stiffness constants.

Non-Parabolic Single Band k · p Model
This model is a simplified version of the multiband model where only a single band is concerned. If one ignores the interaction between the conduction band and the valance bands, which is the case for most III-V direct gap semiconductors, then the conduction band can be described using the following single band parabolic Hamiltonian model [25]: where m * e is the electron effective mass along the growth direction and m * t the transverse one. E c is the unstrained conduction band edge. The fact of using two different effective masses is to account for the anisotropy of the energy band. To increase the accuracy of this model Hamiltonian, the effect of non-parabolicity has to be accounted for, which includes the effect of remote bands in an indirect way without the need to solve the full multi-band Hamiltonian. To this end, an energy-dependent effective mass model has been proposed [26]: where E g is an effective energy gap and does not necessarily match with real energy gap. With this energy dependence, and under the effective mass approximation [27], the eigenvalue problem corresponding to the Hamiltonian in Equation (6) is quadratic and a nonlinear one in terms of energy [28]. To solve such a problem, an iterative algorithm should be invoked to avoid dealing with a fourth-order derivative when ensuring the continuity of probability current across heterointerfaces, which involves only the first derivative of the envelop function.

Finite Element Discretization of a K · P Hamiltonian
The Hamiltonian given above is for bulk semiconductors, and all its terms are position independent such as E c and E v . In addition, the wave vector components k x , k y , and k z are input parameters and no space discretization is needed to calculate bulk dispersion. For bulk semiconductors, the overall band structure is determined by solving the following Schrodinger equation with a Hamiltonian using Equation (1) or Equation (6) in a compact form as: where |ϕ is a column vector for the basis set given in Equation (3) and k is the wave vector. This gives the energy E versus k in different directions, namely k x , k y , and k z . With heterostructures, the band edges become position dependent, i.e., E c (r) and E c (r), and thus constitute the potential profile V(r) of the heterostructure. Furthermore, the Hamiltonian needs an adaption using the envelope function approximation EFA, which is the effective mass theory applied to heterostructures [29].
To show the result of using the EFA as well as the FEM, it is convenient to decompose the bulk Hamiltonian in Equation (1) or Equation (6) into six terms, each one being the coefficient of k 2 x , k 2 y , k 2 z , k x k y , k x k z , and k y k z terms: where the indices run over the coordinates x, y, and z. Then, if we consider the z-axis to be the quantization axis, the adaption of the Hamiltonian includes replacing the wave vector component k z by the corresponding operator −i ∂ ∂z where i = √ −1 , thus the Hamiltonian takes the form: This, as shown shortly, turns the problem of solving the Schrödinger Equation (8) for band structure into a set of simultaneous differential equations for the eigen envelope functions and eigen energies. This formulation can be further extended to any size of Hamiltonian. In our case study, the size is 8 × 8. To preserve hermiticity of the Hamiltonian, the resulting derivative operators in Equation (10) are ordered as follows: This reordering is known as the symmetrization rules [30] given by: which ensure hermiticity of the Hamiltonian as well as the flux continuity across the heterointerface and is equivalent to the Ben-Daniel-Duke boundary conditions in the case of single band model. Within the context of EFA, the wave functions are approximated in terms of envelope functions F µ and basis states |µ : where k t = k x e x + k x e y and r t = xe x + ye y . S is a normalization area of the heterostructures. Thus, the final result of EFA is to solve the following multiband effective mass equation in which the Hamiltonian in Equation (11) acts on an eight-component envelop function vector called spinor F: Here, I is the identity matrix. E k t are the eigen energies of the multiband effective mass equation at in-plane vector k t . It is worth mentioning that, to reduce the size of the problem to be solved, the 8 × 8 Hamiltonian can be block diagonalized into two 4 × 4 Hamiltonians by choosing an appropriate transformation matrix based on a new set of basis states. More details can be found in [31]. To use FEM analysis to solve the multiband effective mass in Equation (15), it is instructive to use either the Galerkin or variational approach, as both lead to the same final expression. In this work, the variational approach is adopted and briefly discussed, and more details can be found in [20]. The variational approach is based on seeking a variational expression corresponding to the set of differential equations to be solved. Following the steps in [20], we end up with the following variational expression: The next step is to minimize this functional using FEM to get the solution of the original problem in Equation (15). To achieve this, the envelop functions are approximated over each finite element as F = N(z) · F = N i F i , where N is an interpolation shape function and F is a nodal column vector with elements representing the values of the envelope function at the element nodal mesh points [16,19]. Substituting F into Equation (15), we obtain the following matrix form to be minimized: which can be cast into a compact form as: where H is a matrix arising from the first and the second terms, and D is a matrix from the last term. Both matrices are of size 8 · n s × 8 · n s where n s is the number of degrees of freedom used in the system.
The minimization procedure with respect to the unknown nodal values of F requires ∂V exp ∂F † = 0, which yields: This is a generalized eigen value problem for the energies E and the representative nodal values of the envelope functions at the mesh points. In constructing the global matrix H and D, an element matrix overlay approach is usually used, which requires considerable care since we are dealing with up to eight simultaneous equations. The problem in Equation (19) can be efficiently solved using well-established algorithms devoted for large eigenvalue problem.

Simulation Results and Discussion
In this section, three applications are shown. First, a case study was solved to validate the accuracy of the simulations to published results using an energy-dependent effective mass non-parabolic model on a structure of a single quantum well [28]. This is a 100 Å GaAs thin well layer surrounded by AlGaAs barrier layers. The barrier height was set to 0.276 eV and the effective masses of electrons at the band edges were taken to be m w /m 0 = 0.067 in the well and m b /m 0 = 0.09 in the barrier. Here, the subscripts w and b denote the well and barrier, respectively, and m 0 is the electron free mass. The in-plane wavevector components k x and k y were set to zero for the purpose of comparison. Using this model, Equations (6) and (7) under the EFA were solved repeatedly by updating the effective masses and eigen energies until their convergence. The results of the second-order Lagrange finite element method are shown in Table 1 along with the results of two other numerical methods, namely the Finite Difference Method (FDM) [32] and the Transfer Matrix Method (TMM) [33], which are also included in the simulation tools developed. The calculation was carried for both the parabolic and non-parabolic case and the results are compared with Samir et al. [28]. Table 1 lists the confined eigen energies in quantum wells of different widths ranging from 5 Å to 200 Å. The results are in excellent agreement with the reported data. The results also clearly show that the single-band parabolic model leads to energies slightly over estimated for higher lying states. Therefore, our results confirm the effects of including non-parabolicity, which shrinks the confined energy gaps with a noticeable shift in the higher energies.
An illustrative example for such effect is given in Figure 1, where the probability density of the first five states of a 100 Å GaAs quantum well embedded in AlGaAs is shown for the parabolic and non-parabolic case. In both panels, the energies E n are indicated by the thin horizontal lines, and the probability densities |ϕ(z)| 2 by the shaded colored areas. There are only three confined states for electrons. The tunneling of the eigenfunctions through the barriers is larger for states with high energies and decreases from the parabolic to non-parabolic. The second application that is investigated in this work deals with the calculation of the bulk band structure of InAs material strained by a GaAs substrate of 5.6533 Å lattice constant. Due to the lattice mismatch between these two material, the grown bulk layer experiences an intrinsic compressive strain of −7.09%. The value of strain used in this case is calculated using the formula (a sub − a)/a × 100% where a sub and a are the lattice constant of the substrate and of InAs material. The overall conduction band and valence bands dispersion using the 8 × 8 k · p Hamiltonian in Equation (1) is obtained by diagonalizing Equation (1) or equivalently by solving the following secular equation at each wavevector k ≡ k x , k y , k z : In Figure 2, the overall band structure curves of strained bulk InAs and unstrained standalone bulk InAs are shown along the k x and k z directions, respectively. These curves are dispersions for the conduction band (CB), Heavy Hole (HH), Light Hole (LH), and spin-orbit split-off (SO) valance bands. In Figure 2 (right), the 6 × 6 k · p Hamiltonian curves are shown. This 6 × 6 Hamiltonian couples only the three valance bands and ignores the coupling with conduction band, which is justified only in wide energy gap materials. Within this 6 × 6 Hamiltonian, the conduction-band states, which are not shown in this figure, can be treated separately using single-band model. In narrow energy gap materials such as InAs, the 6 × 6 Hamiltonian results show a discrepancy for higher in-plane wave vectors. Figure 3 shows the band structures of unstrained bulk InA along crystallographic directions pointing to high symmetry points. The Hamiltonian used is the 14 × 14 one from [34], which extends the interaction to more higher bands. The results confirm that including the conduction band into the Hamiltonian increases the conduction band for higher k x values while pushing down LH and SO bands. This effect is negligible on HH band. In Figure 3 (left), the conduction band is parabolic in both directions k x and k z while the valence bands exhibit a non-parabolicity due to mixing effects. The strain lifts the HH and LH bands degeneracy at the zone center. The parabolic isotropy is in fact a reflection of the s-like character of conduction-band states in a bulk semiconductor. Table 2 lists the k · p band parameters used for the calculation [7,35]. It should be noted that the k · p method results are sensitive to the input parameters.
Moreover, the k · p solutions are subject to a well-known problem of unphysical "spurious" solutions, especially with Hamiltonians that include more remote bands or when band parameters have large difference across layers' interfaces [31].   Such unphysical solutions should be identified and removed from the region of interest. The spurious solutions arise from the incompleteness of the basis states used for the construction of the Hamiltonian [37]. This is still a standing problem; different approaches and models have been proposed to eliminate them but each approach has its drawbacks. Such solutions may appear as highly oscillatory or strongly localized wave functions within the band gap or dispersion curves bending in the wrong direction [38][39][40]. Three main reasons have been identified for the source of the spurious solutions. First, the input experimental parameters can trigger wing-band solutions within the band gap originating from second-order Hamiltonian's terms at large-k values [38,40,41]. These solutions can be avoided by parameter rescaling procedures [42,43]. Second, the interface between two materials can be subject to localized spurious solutions due to improper boundary condition for the envelope function components. These solutions can be removed by a proper ordering of differential operators [43,44]. Finally, the spurious solutions may also be triggered by the chosen discretization scheme and the mesh width, and it was found to be related to the ill-representation of the first-order derivative terms of the Hamiltonian [40,45].
In the developed simulation tools, two elimination models are implemented. The first one is based on adding a small correction to the off-diagonal elements of the Hamiltonian [41] to discard terms responsible for the spurious solutions and the second one modifies Kane's parameter P cv [37,38].
Finally, the last test case illustrates the ability of the tools to model a realistic optoelectronic device where the miniband dispersions of the active region of an infrared photodetectors based on type-II InAs/GaSb superlattices are calculated. The studied superlattice is mainly composed of nine InAs/GaSb periods with thicknesses around 44 Å/21 Å, respectively. The InAs layers are chosen to be thicker than GaSb layers to efficiently reduce the dark current [46]. This superlattice is designed to operate in a long wavelength infrared radiation (LWIR) 8-14 µm atmospheric detection window. The potential profile V(r) of this active region is periodic and this periodicity is handled using the Bloch theorem, which greatly simplifies the computational task. Within this theorem, the wave function in two adjacent superlattice periods differs only by a phase shift: ϕ(z) = ϕ(z − p)e ik p ·p where k p is the Bloch wavenumber and p is the superlattice period. From a computational point of view, this boundary condition is implemented by overlaying the last element matrix onto the first one modulated by the factors e ±ik p ·p in the global matrix system. The Dirichlet boundary condition is implemented by setting the envelope function to zero at the external borders of the structure. Type II InAs/GaSb superlattice is made of narrow energy gap materials in which the coupling between conduction bands and valance bands is strong; therefore, the 8 × 8 k · p Hamiltonian in Equation (1) should be used. To account for this strong interaction, the original electron effective mass m * e in the diagonal Hamiltonian matrix as well as the original Luttinger parameters γ 1 , γ 2 , and γ 3 should be updated according to: where m * c , γ L 1 , γ L 2 , and γ L 3 are now the original effective mass and the original Luttinger parameters, respectively. These normalizations are only required with Hamiltonians of size greater than 8 × 8. This is because the strong interaction is explicitly included in the formulation of the eight-band model, whereas, in smaller size k · p models, the valence and conduction bands are decoupled and the Luttinger parameters are adjusted using Löwdin perturbation method [47]. Using the Löwdin perturbation, method, the Hamiltonian is reduced to a smaller number of relevant bands, including the effect of the remote bands perturbatively. In all calculations, the symmetrization rules given by Equation (12) are applied to taken into account the real space dependence of the bulk band parameters.
The band edge profile potential emanating from a structure of a nine-period InAs/GaSb44 Å/21 Å superlattice, simulated using our model is shown in Figure 4 (right). The corresponding miniband energy-dispersions versus in-plane wavevector k t solved using Dirichlet boundary conditions are shown in Figure 4 (left). In this figure, the reference of energy is taken to be the unstrained conduction band edge of InAs layers with a VBO value of 560 meV. The formation of minibands confirms the quantum mechanical intuition that, since electrons and holes in InAs/GaSb Type II superlattice are separated in space, confined electrons in thin InAs layers form their mini-bands due to overlap of electron wave functions. Similarly, the interaction among holes in thin GaSb layers leads to hole minibands. The Hamiltonian used is the full 8 × 8 k · p model without block diagonalization. The fundamental electron and heavy-hole envelope-function component moduli squared for nine periods solved using Dirichlet boundary conditions are shown in Figure 4 (right). At k t = 0, a classification of minibands states as heavy-hole, light-hole, and split-off states is possible, and are denoted by HHn, LHn, and SOn (n is the quantum number of the subband), respectively. The minibands nomination are assigned based on the dominance of envelope functions at the band edges [7], although this is usually not achievable with the eight band model and with all other higher multi-band k · p models exhibiting strong band mixing. The case is worse with narrow gap systems such as InAs/GaSb and type II heterostructures with a broken gap. The distribution of the eight components of the envelope function spinor F µ |F µ for the 21 eigenstates solution of Equation (19) near the effective gap is plotted in Figure 5. In this figure, the states numbered 4 and 20 have the largest components F µ |F µ among all other states. In addition, State 20 has the largest component F 1 |F 1 corresponding to S|S , indicating a s-like character of conduction-band, while State 4 is mostly dominated by the two largest components F 2 |F 2 and F 6 |F 6 , indicating a mixture of valence band character. Therefore, these two states correspond, respectively, to the lowest state of a conduction miniband C1 and a highest state of a valence miniband HH1. The band dispersion of the SOn states are of less practical interest due to their lower energy range compared to the range of HHn and LHn states and therefore are not shown. The energy dispersion relations and the corresponding fundamental electron and heavy-hole envelope-function component moduli squared for a nine-period InAs/GaSb 44 Å/21 Å grown on GaSb substrate solved using Periodic boundary conditions are illustrated in Figure 6. Since there are nine degenerate set of subbands under Dirichlet boundary conditions, each of which belongs to one period, the strong interaction lifts the degeneracy into nine states for each subband, as shown in Figure 4. In general, each state splits into N p states where N p is the number of periods. Using periodic boundary conditions, because of the double degeneracy of ±k p states, E(+k p ) = E(−k p ), the degeneracy is lifted only into five states for each subband, as shown in Figure 6. For both cases of boundary conditions, the dispersion relations converge at higher in-plane wavevector k t as the interaction becomes weaker. The time for the entire calculation in this case study using coarse meshing is 608.18 s. Obviously, the higher is the complexity of design, the longer is the execution time. This time increases with the increase in the number of layers, the size of the Hamiltonian, and the number of in-plane wave vector sampling points. However, it could be optimized with the FEM approach where, for larger layers, one can use coarse mesh instead of fine mesh for thin layers. In fact, it could be further optimized by using Hermit shape functions.
The full form of the 8 × 8 Hamiltonian used in this work is reflected in the anisotropy of band structure. Figure 7 shows the fundamental heavy hole sub-band structure in polar coordinates for both full and block diagonalized Hamiltonian. The anisotropy is strong in the transversal [110] direction for the full form and is absent in the block diagonalized one, which yields unrealistic isotropic in-plane masses in the later case. This anisotropy is usually neglected in previous works where the axial approximation is usually used after block-diagonalization. The transversal directions along which the dispersion relations exhibit the largest difference are [100] and [110]. This anisotropy of the sub-band is estimated by the difference E [100] (k t ) − E [110] (k t ). Between (k x = 1, k y = 0) (Å −1 ) and (k x = 1, k y = 1) (Å −1 ), it was found to be about −90 meV, which indicates that the effective mass is slightly larger along the [110] direction than the [100] direction. This observation confirms that the hole sub-bands are neither parabolic nor isotropic, as can be seen from the dispersion curves.  Despite the difference between the two types of boundary conditions used, the effective band gaps and band edges at the zone center are expected to be comparable since there is no phase shift between the envelope functions at k p = 0. As indicated in Figures 4 and 6, the effective bandgap is determined by the difference between the bottom of the lowest electron miniband (C1) and the top of the highest hole miniband (HH1) and it has the value of 0.13 eV, which is in excellent agreement with the one experimentally reported by [48] and theoretically with the absorption spectrum reported in Figure 8. This effective gap can be tuned by changing the layer thickness and is primarily determined by the position of the C1 bottom edge since the HH1 band width is less sensitive to layer thickness [49]. The HH1 band width is too small which might be interpreted as it experiences a balanced effect between higher and lower bands. The absorption spectrum at photon energyhw is derived from Fermi's golden rule as: where ε 0 and C are the speed of light and permittivity in a vacuum, and n r and L are the refractive index and the thickness of the superlattice, respectively. f n c and f m v are the Fermi functions for the mth valence sub-band with an energy E m v and the nth conduction sub-band with an energy E n c , respectively, given by: To account for different scattering relaxations causing finite transition linewidth, a normalized Gaussian or Lorentzian L(k t ,hw) distribution function with a broadening γ parameter can be chosen. The summation over k t takes into account the spin degeneracy. The calculation of the optical transition matrix elements Ψ nk t |ê · P|Ψ mk t between the hole sub-bands and the electron subbands for TE (e⊥ẑ)or TM (e ẑ)-polarized optical beams requires the evaluation of a dense 8 × 8 momentum matrix P: Then, the final expression for the momentum matrix elements reads: To accelerate the evaluation of these elements, the derivatives of the Hamiltonian elements are calculated explicitly in terms of spinor envelope components and are given in Appendix A. Using Equations (A1) and (A2), the average time required to calculate the matrix element of a single transition at a given in-plane wavevector is reduced from 7 s to 0.7 s compared to the direct Hamiltonian derivation method. Figure 8 shows the total absorption spectrum for nine-period InAs/GaSb 44 Å/21 Å type II superlattice using both kind of boundary conditions. The cutoff energies for different types of interband transitions are consistent with the effective band gaps between different sub-bands predicted in Figures 4 and 6. The maxima observed in Figure 8 of around 550 meV indicate a large wavefunction overlap of around this energy value due to larger in-plane wavevector states' contributions, where the hole dispersions are neither parabolic nor isotropic. Figure 9 shows the largest envelope function spinor component dependent on in-plane wavevectors near the LH1 miniband. Clearly, the in-plane states' contributions are high to produce the absorption maxima. To clarify the sensitivity to the input parameters in Figure 10, the absorption spectrum for different valence band offset (VBO) values is calculated. As could be inferred from Figure 10, an increase of InAs/GaSb superlattice's VBO from 510 meV to 575 meV leads to an increase of the absorption coefficient. Such a sensitivity indicates that the type II superlattice structure should be properly designed to achieve maximum quantum efficiency in real devices such as infrared detectors. In addition, for the LWIR window, a shift of nearly 1 µm in the cutoff wavelength is due to 5 meV change in the VBO. Our calculations show that the LWIR cutoff wavelength is in a good agreement with the experimental data from reference [48]. There exists a further complication with InAs/GaSb superlattices. Although the separation of the confined states and minibands can be tuned to a desired IR windows by adjusting the InAs and GaSb layer thicknesses, it was found that the physical and optical properties and device performance depend critically on the composition and structure of interfaces formed at the heterojunctions. Recent studies demonstrated that the interface InAs/GaSb is not always abrupt and in fact has chemical asymmetries in compositional profile due to segregation of Sb and In atoms at the interfaces [50,51]. This segregation leads to the formation of unintentional interfaces (IFs) with different composition configurations. For InAs/GaSb SLs, IFs often turn out to be InSb-like or GaAs-like. Previous studies modeled IF's effect either by fine-tuning the VBO values or the use of an interface potential profile (short-range delta function potentials centered at the interface [52] and graded or asymmetric interface [53,54]). However, these approaches miss the effect of the intrinsic strained property of IFs on the electronic structure and the optical properties. The most successful approach is to insert InSb layers to account for large IFs lattice mismatch and the strain effect. A structure with an IF layer is obtained by introducing an interface InSb layer between InAs and GaSb on both sides, while the thickness of one superlattice period is kept fixed. The change in the fundamental effective gap as a function of the number of periods N are plotted in Figure 11 for InAs/GaSb 44 Å/21 Å SL without and with InSb IF layer forced at 1.2 Å.
For small number of periods, the effective gap exhibits a strong variation, which indicates the quantum mechanical interaction and the start of miniband formation, while, for a higher number of periods, the effective gap tends to converge to the fundamental minibands gap of the superlattice. For very small periods, the large scatter of eigen energies is a signature of the single quantum well confinement. The presence of InSb IF layer induces a decrease in the effective energy about 25%. This is because the valence band state energies are shifted upward due to the presence of hole states in InSb layer lying in a higher level than those in InAs and GaSb, whereas this effect is negligible with conduction band states since these states in InSb are in between the corresponding ones in InAs and GaSb. The overall effect of including IF layers on the HH miniband is minor with an upward shift less than 6 meV. The absorption coefficients with effect of interfacial layers taken into account are shown in Figure 12 indicating an increasing shift in the cut-off wavelengths with IF layers thickness. Figure 11. Variation of the inter-miniband band transition energies at k x = k y = 0 with the number of periods for InAs/GaSb 44 Å/21 Å type II superlattice. The HH1-C1 and LH1-HH1 transitions are denoted by brown and mauve shaded zones, respectively. The effect of including InSb interfacial layers are also shown (left). The blue square markers are energies of states above the bottom of C1 miniband edge and below the top of LH1 miniband edge, while those at middle are for HH1 miniband. The value used for VBO is 560 meV.

Conclusions
In this work, an efficient, accurate and rigorous model for the overall band structure calculation based on the k · p approach is presented and tested for cubic zincblende bulk semiconductor materials as well as for their heterostructures. These modeling tools are useful for the design of optoelectronic devices and the estimation of their performance. The tools include different k · p Hamiltonian formulations namely: 14 × 14, 8 × 8, 6 × 6, 4 × 4, and single band non-parabolic model. These models take into account the strain effects and can be invoked easily within the modeling steps. The results of the using a non-parabolic model are in excellent agreement with published data. The predicted cutoff energies for different types of interband transitions are consistent with the calculated band dispersion curves of InAs/GaSb Type II broken gap superlattice and demonstrate excellent agreement with the experimental data reported in [48]. Our modeling shows that the IF's lattice mismatch and induced strain yields a significant increase in the absorption coefficient as well as in cut-off wavelengths. As a result, a wide tunable range of optical properties can be obtained with IF layers. Finally, the model can be extended easily to cover materials belonging to wurtzite crystal symmetry just by generating the correct Hamiltonian.
Author Contributions: Z.B. and A.B. were mainly responsible for writing the original draft, developing and validating the model and analysing the important information presented in the manuscript. N.S. contributed to the funding acquisition and project administration.

Acknowledgments:
The authors acknowledge the financial support of the Ministry of Higher Education and Scientific Research, Algeria. The authors are also grateful to the valuable discussion with Phil Picton that contributed to enriching the final version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Momentum Matrix Elements
We provide here the explicit derivation of matrix elements nPm = Ψ nk t |(x +ŷ) · P|Ψ mk t corresponding to TE transition mode between the initial spinor |Ψ mk t and the final spinor |Ψ nk t . These are for the 8 × 8 KP Hamiltonian given by Equation (1) and are evaluated for any direction in transversal reciprocal k t (k x , k y ) space and in terms of spinor envelope components |S mk t from the initial vector state and |F nk t from the final vector state. In the following, the subscript k t has been dropped for convenience. (A2)