1. Introduction

Semiconductor superlattices (SLs) are semiconductor materials that incorporate a one-dimensional periodic potential with a period larger than the crystal lattice^[1]. In these materials, the motion of charge carriers is confined along one of the dimensions, resulting in unique quantum properties and modification of the band structure of the host semiconductors. Even several decades after the first study of SL was published, scientific groups still discover new applications^{[2, 3]} and curious properties^{[4, 5]} of these quantum structures.

One of their most promising applications is the creation of THz-emitting devices. The THz radiation has a variety of essential applications in medicine, astronomy, and telecommunications^{[6, 7]}, which indicates the need to produce new cutting-edge tunable THz sources operating at room temperature. Multiperiod SLs are the prime candidates for this job because they have all the components required for light amplification and their designs seems to be simpler and more reliable, than designs of quantum cascade lasers (QCLs)^{[8, 9]}. The ongoing simplification of designs of QCLs towards the SLs^{[10, 11]} confirms this conclusion.

Experimental studies and fabrication of SLs have become possible in the last several decades due to significant progress in growth technologies. The most valuable contribution to this research field was made by molecular beam epitaxy (MBE), which brought to life more and more devices with ultra-thin layers. Moreover, a new type of SL, super-multiperiod (SMP) structures, i.e., SLs with a number of layers from several hundred up to few thousands have been recently implemented by MBE^{[9, 12−16]}.

In both cases (ordinary or SMP SLs), researchers need to solve the Schrödinger equation with the given periodic electronic potential. The introduction of such potential leads to the transformation of the initial electronic Bloch states and minibands formation. Currently, there are plenty of numerical methods to get this information qualitatively and quantitatively. But they all incorporate a variety of assumptions, limitations, and approximations. For instance, kp-methods consider states in the vicinity of the selected point of the Brillouin zone, whereas the density functional theory calculations require a complex and intense computation procedure^[17]. At the same time, all these methods are the most appropriate ways to tackle the problem of determination of the electronic properties of SL. However, in more intriguing cases, the theory is yet to be solidified, e.g., SMP SLs under external electric fields. One of the vivid examples of this fact is the problem of predicted low gain values in SLs^[11]. It was suggested, that in similar SL structures, the gain is strongly affected by elastic scattering for a higher temperature. To the best of our knowledge, such comprehensive investigations, which considers combination of both experimental and theoretical approaches, have not been conducted yet.

In this work, a new theoretical approach for the modelling of SMP SL parameters is developed. Our method combines theoretical and numerical models, including 8 band kp- and shooting methods, as well as Monte Carlo (MC) simulations. The developed approach was proved by the study of the properties of Al_xGa_1-xAs/GaAs SMP SL samples of the finest quality obtained by MBE. These samples were also examined by the express method of photoreflectance spectroscopy (PRS) and high-resolution transmission electron microscopy (HR ТЕМ). Conducted research establishes the accuracy, applicability, and versatility of the developed approach.

2. Experimental methods and structures

The Al_xGa_1-xAs/GaAs SMP SLs of high quality were grown on one-sided semi-insulating GaAs (100) substrates using the Riber 49 MBE setup. The reflection high-energy electron diffraction system was used to control the growth processes in situ. The thickness of the GaAs layers in the synthesized SMP SL structures differs across the samples, while the barrier Al_xGa_1-xAs thickness remains constant. The Al mole fraction x in the barrier layers was equal to 30%. The grown SMP SL structures were doped identically in barriers and wells with silicon as an n-type dopant. All samples had the same concentrations of n-type donor dopant (n = 10¹⁶ cm⁻³). Table 1 contains all the essential parameters for the samples.

Obtained samples exhibit the excellent quality of the SMP SLs which was confirmed by different methods. The detailed metrology of the grown doped and undoped SMP SL Al_xGa_1-xAs/GaAs structures with period numbers of 100–400 was carried out by high-resolution X-ray reflectometric (HRXRR), deep X-ray reflectometric (DXRR) and high-resolution diffractometric (HRXRD) studies, in particular, using a synchrotron radiation source^{[12, 13, 15, 16]}. The difference in the thicknesses of the grown layers, estimated by three methods and over the entire depth of the optimized SL structures, was no more than 0.5–2.0 percent, the interface root mean square (RMS) roughness was a few Å, and the composition inaccuracy was ~0.5 atomic percent. The high sharpness of SLs’ interfaces was also demonstrated by high-resolution transmission electron microscopy (HR TEM) (see Fig. A1 and Supplementary materials)^{[13, 18]}.

One of the main methods used for precise energy spectrum retrieving was PRS^{[19, 20]}. PRS is the non-contact, non-destructive, and highly sensitive express method for diagnosing nanostructures. It allowed us to determine the energies of interband optical transitions in nanostructures not only at low^[21], but also at room temperature ranges^{[22, 23]}. One can also obtain information about higher-order optical transitions in different nanostructures^{[24, 25]}, including type-Ⅱ SLs^[26]. The high quality of the grown samples makes this method applicable for the precise examination of optical properties of SMP SLs even at room temperatures. PRS is based on the application of both the probe laser beam, which varies in a wide energy range, and the periodic modulated laser beam. The electron−hole pairs photoinduced as a result of laser exposure are separated by the near-surface built-in electric field. At the frequency of the periodic laser perturbation, a change in the reflectance coefficient is synchronously detected, which underlies the recorded signal^[19]. During the study of low-dimensional semiconductor structures, the formation of a signal can occur due to a change in the position of quantum-dimensional levels under the action of a change in the built-in electric field due to the Stark effect^[27]. In case of SLs, the quantum-dimensional levels can be transformed into several energy minibands of finite width. Therefore, singularities in the photoreflectance (PR) spectra can be observed at energies corresponding to the maximum and minimum energies of possible transitions between minibands (near the Г and П points of the Brillouin zone of the SL)^[28−30].

The PR spectra were measured using an experimental setup assembled on the basis of an IKS-31 diffraction spectrometer. The source of modulating radiation was a highly stable SSP-DHS-405 violet laser with an operating wavelength of λ = 405 nm. The laser beam was aimed at the exact area of the sample where the spectrometer's monochromatic light was directed to. Modulation signal was recorded using a silicon photodiode, a preamplifier, and a lock-in amplifier (SR830, Stanford Research). The measured PR spectra represent a series of differential signals, the inflection points of which determine the energy values of optical transitions^[19]. To get accurate energy values and characterize spectral features, PR spectra were transformed using the Kramers–Kronig relations (Kramers–Kronig method, KKM)^[31]. The Appendix B provides a detailed explanation of this transformation. The extrema of transformed signals determine the precise energy values of the observed optical transitions.

3. Results and discussion

The presented SMP SL structures in Table 1 were initially designed for the light amplifiers’ implementation. The geometry of these SLs was inspired by the structures proposed in Ref. [8] and underwent a few modifications for output characteristics improvement, including estimated gain^[9]. When the electric field is applied to such SLs, light amplification becomes possible due to the emergence of optical intersubband transitions. Since these transitions occur between energy minibands of SL, their relative position, i.e., the emitting wavelength, can be tuned with the applied field. In this case, all output characteristics of the SLs depend on their parameters (layer thicknesses, composition, and doping concentrations). For instance, thinner Al_0.3Ga_0.7As barriers maintain higher dipole matrix element values for optical transitions, resulting in higher gain values. Such qualitative considerations provide a general idea of the structure design but do not define it completely. Therefore, to fully define the structure, i.e., to choose structure parameters suitable for the task, it is necessary to carry out numerical simulations of the SL energy spectrum.

3.1 Initial design and simulation approaches

The initial designing process was based entirely on the results of the classical shooting method (CSM) simulations^[32]. Applied CSM utilizes the approximation of the single-band kp-model, an easily-implementable computation scheme which showed a good agreement with experimental results for a single quantum well, multiple quantum well structures, and ordinary SLs^{[33, 34−37]}. This method has been also selected because it is frequently utilized to address similar QCL simulation tasks^{[10, 38]}. CSM retrieves the electronic eigenstates, which define the electronic energy spectrum for a given SL periodic potential. A boundary value problem for the Schrödinger equation in CSM is reduced to the Cauchy problem, i.e., an initial value problem^[32]. The electronic eigenstates are determined by 'shooting' solutions from one boundary condition with different initial conditions until the other boundary condition is met. Once the eigenvalues are found, one can employ a transfer matrix approach to retrieve wavefunctions via the propagation procedure. The transition energies are defined as the difference between the values of the corresponding eigenenergies (minibands' boundaries). The electronic energetic spectrum is then utilized in a stochastic MC approach to gather information about the kinetics of the SLs, such as transition rates and gain estimations. The Fermi golden rule is applied to calculate the scattering rates used in MC, including phonon and impurity scattering rates. For example, to calculate the impurity scattering rate the following formula can be employed:

$${W}_{ik,j}=\frac{{m}_{j}^{\left|\right|}{e}^{4}}{4\pi {\epsilon }^{2}{\hslash }^{3}}\underset{0}{\overset{\pi }{\int }}\frac{{F}_{ij}\left(q\right)\rm{d}\theta }{{q}^{2}},$$

(1)

where $q(\theta ) = |k-{k}'|,{F}_{ij} = \int \rm{d}{z}'{n}_\mathrm{D}({z}'){[\int {\psi }_{i}(z){\psi }_{j}^{*}(z)\mathrm{exp}(-q|z - {z}'|)dz]}^{2}$, $k,{k}'$ are plane wave vectors, $\theta$ is the angle between $k,\ k'$, $V$ is impurity potential, ${n}_\mathrm{D}$ is doping sheet density for the considered layer, $i,\ j$ are indexes of eigenenergy levels of SL between which the scattering occurs, ${m}_{j}^{\left|\right|}$ is in-plane electron effective mass for the j-th energy level, $e$ is the elementary charge, the dielectric permittivity $\epsilon ={\epsilon }_{0}{\epsilon }_{r, 0}$, ${\epsilon }_{0}$ is the vacuum permittivity, ${\epsilon }_{r, 0}$ is the static dielectric constant. For scattering on longitudinal optical (LO) phonon:

$${W}_{ik,j}=\dfrac{{m}_{j}^{\left|\right|}{e}^{2}{\omega }_\mathrm{LO}({\epsilon }_{r,\infty }^{-1}-{\epsilon }_{r,0}^{-1})\left({N}_\mathrm{ph}+\dfrac{1}{2}+\dfrac{1}{2}\right)}{8{\pi }^{2}{\hslash }^{2}{\epsilon }_{0}}\underset{0}{\overset{2\pi }{\int }}{F}_{ij}\rm{d}\theta ,$$

(2)

where $q=\left|k-k'\right|,{F}_{ij}=\dfrac{\pi }{q}\int \rm{d}z\rm{d}{z}^{'}{n}_\mathrm{D}\left({z}^{'}\right){\psi }_{j}\left(z\right){\psi }_{j}\left({z}^{'}\right){\psi }_{i}\left({z}\right){\psi }_{i}^{*}\left({z}^{'}\right)\cdot \mathrm{exp}\left(-q\left|z-{z}^{'}\right|\right)$$, \omega_{\mathrm{LO}}$ is the LO phonon frequency at $Q=0,Q$ is the phonon wave vector, $N_{\mathrm{ph}}$ is the mode independent phonon occupation number, ${\epsilon }_{r, \infty }$ is the high frequency dielectric constant.

To estimate gain based on dipole matrix element $D$ and doping concentrating ${n}_\rm{e}$, the following formula, obtained from two level model was used^[32]:

$${\alpha }_\rm{g}=\frac{4\pi \omega {e}^{2}{n}_\rm{e}}{c{n}_{\omega }}\frac{\gamma {D}^{2}}{{\left({E}_{21}-\omega \right)}^{2}+{\gamma }^{2}}\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(\frac{{E}_{21}}{2kT}\right) ,$$

(3)

where $\gamma ={\tau }_{1}^{-1}+{\tau }_{2}^{-1}$, $D$ is the dipole matrix element, ${E}_{21}$ is the energy of transition, $\omega$ is the angular frequency of radiation, ${n}_{\omega }$ is the refractive index of the structure, $T$ is the electron temperature, $c$ is the vacuum speed of light, ${\tau }_{i}$ is the lifetime of the level $i$. The influence of the coherence of the electron transport^[39] was considered implicitly in the MC approach via collisional broadening in the $\gamma$ parameter^[40]. Usually, a depolarization shift is also taken into account in classical designs of THz QCLs^[41]. However, this shift was neglected in that case due to several reasons. First, in case of SMP SLs with the considered doping levels and slightly different lasing scheme (compared to QCLs), further investigation is still required to analyze the impact of such depolarization shift. Moreover, all the computations were performed for room temperature and biases of approximately 50 kV/cm. There is experimental evidence that in such cases, the depolarization effect influences the shape of the gain curve, but the shift becomes less pronounced^{[42, 43]}. Also, the figure of merit $F_{\mathrm{M}}$ was utilized to take into account both dipole matrix element and population inversion between upper $u$ and lower $l$ states^[44]:

$${F}_\mathrm{M}={D}_\rm{ul}^{2}{\tau }_\rm{eff}={D}_\rm{ul}{\tau }_\rm{u}\left(1-\frac{{\tau }_\rm{l}}{{\tau }_\rm{ul}}\right),$$

(4)

where ${D}_\rm{ul}$ is the dipole matrix element for transition between upper and lower levels, ${\tau }_\rm{ul}$ is the inverse transition rate between upper and lower levels.

During the design phase, we varied the thicknesses of all the layers and layer composition to maximize both the figure of merit and gain estimation while preserving the emitting wavelength in the THz range. In order to check the calculated design, a series of SMP SL Al_0.3Ga_0.7As/GaAs structures were grown (Table 1). SMP SL structures had sharp interfaces, their thicknesses and composition deviations were minimal (see Fig. A1 and Supplementary materials). The creation of such high-quality SMP SLs with required layer thicknesses and acceptable interface widths became possible due to laboratory and synchrotron studies^[12], as it was mentioned in the previous section.

The optical properties of the samples were examined by PRS to determine the electronic energy spectrum. The PR examination is a suitable verification tool because the electronic energy spectrum and the optical output features of the structures are closely interconnected. Therefore, for all the grown SMP SL structures, PR spectra were measured and transformed according to the method described in the previous section and Supplementary materials. The data obtained were compared with CSM results. For instance, the e1−hh1 transition energy, i.e., the difference between the first electronic (e1) and the first heavy hole (hh1) energy levels, was compared with the position of the first peak in the transformed spectrum.

We anticipated only minor deviations between experimental and simulated values since CSM is widely applied in similar QCL design tasks^[32]. However, the PR analysis (Fig. 1) of sample no. 1 showed a significant deviation from the expected band structure: simulated energy values (green arrows in Fig. 1) and experimental values (black dotted arrows in Fig. 1) differed for e1−hh1 and e1−lh1 transitions by 15 meV and 23 meV, respectively. That results in a total error of about 25% and makes CSM inapplicable for the proper and complete electronic energy spectra analysis of SMP SLs.

As CSM was not suitable for the analysis of the transformed PR spectrum, we had to find an alternative simulation approach. The first idea was to test more advanced purely numerical methods: the finite elements method (FEM), which numerically solves the Schrödinger equation for an infinite SL on the specified grid^[32]; the transfer matrix approach with the Airy functions and immense modifications of the shooting method, which also provides the numerical procedure to obtain the minibands position of SL^{[32, 37]}. We reduced the error to 10%, but the accuracy of these methods was still insufficient for an unambiguous analysis of experimental data. Therefore, we changed our simulation approach to a more elaborate one.

It is well known that the theoretical study of semiconductor heterostructures and SLs can be done in various ways: the empirical pseudopotential method (EPM)^{[33, 45−48]}, the tight binding (TB) method^{[33, 48]}, and the kp-method^[33−36]. All of the mentioned methods allow to simulate SLs with a large number of periods (here, SMP SLs). However, there are few works only^{[8, 11]} that include the experimental study of such SLs combined with theoretical and numerical analysis. Some of these methods even assume that the SL under consideration is infinite^{[33−36, 46−52]}. This is also an approximation to be verified with the experimental data. Taking into account all of these facts, it was not clear which method to choose. The EPM has an advantage over TB, since EPM requires fewer adjustable parameters and a suitable basis for the solving of the Schrödinger equation can always be chosen. The EPM makes it possible to obtain results those are in good agreement with the experiment^{[33, 45, 46]}, since it can account any number of bands contributing to the formation of energy states. However, such many band EPM considerations impose additional computational requirements even for the structures whose dimensions are just a few nanometers. Although there are ways to get around these difficulties^[33], the computers used have not allowed one to study SMP SLs with a supercell of more than ~2–3 nm in size by this method. The kp-method is suitable for investigation of nanostructures larger than a few nanometers, including SMP SLs, as it was demonstrated in Refs. [36, 50]. This method is also among the most appropriate methods for analyzing experimental data. The complexity of the kp-method increases with the number of bands considered in the calculation^{[34−36, 53]}. The initial simulation approach included only an approximation of the single band kp-method, which probably was the reason for poor accuracy. Thus, we selected the below-described 8-band kp-method, also known as the Kane model, for the appropriate theoretical study of the SMP SLs.

3.2 Theoretical study of SMP SLs: the Kane model

The 8-band Kane model was previously used to calculate the energy spectrum and interband matrix elements of similar heterostructures and SLs^[49]. One of the key advantages of the Kane model is that it allows to consider the nonparabolicity of the charge carriers' dispersion law. Another bonus is the insignificant computational burden of the model, even in our case of SMP SLs^[35]. The model can be formulated as an eigenenergy problem for the Schrödinger equation with the following Hamiltonian:

$$\begin{aligned} {H}_{\mathbf{k}{\mathbf{k}}'}\left(\mathbf{K}\right)=\left(\begin{array}{ccc}{H}_{\mathbf{k}{\mathbf{k}}'}^{\left({\Gamma }_{6}^{c}\right)}\left(\mathbf{K}\right)& {H}_{\mathbf{k}{\mathbf{k}}'}^{{\Gamma }_{6}^{c}{\Gamma }_{8}^{v}}\left(\mathbf{K}\right)& {H}_{\mathbf{k}{\mathbf{k}}'}^{{\Gamma }_{6}^{c}{\Gamma }_{7}^{v}}\left(\mathbf{K}\right)\\ {\left[{H}_{\mathbf{k}{\mathbf{k}}'}^{{\Gamma }_{6}^{c}{\Gamma }_{8}^{v}}\left(\mathbf{K}\right)\right]}^{\dagger}& {H}_{\mathbf{k}{\mathbf{k}}'}^{{(\Gamma }_{8}^{v})}\left(\mathbf{K}\right)& {H}_{\mathbf{k}{\mathbf{k}}'}^{{\Gamma }_{8}^{v}{\Gamma }_{7}^{v}}\left(\mathbf{K}\right)\\ {\left[{H}_{\mathbf{k}{\mathbf{k}}'}^{{\Gamma }_{6}^{c}{\Gamma }_{7}^{v}}\left(\mathbf{K}\right)\right]}^{\dagger}& {\left[{H}_{\mathbf{k}{\mathbf{k}}'}^{{\Gamma }_{8}^{v}{\Gamma }_{7}^{v}}\left(\mathbf{K}\right)\right]}^{\dagger}& {H}_{\mathbf{k}{\mathbf{k}}'}^{\left({\Gamma }_{7}^{v}\right)}\left(\mathbf{K}\right)\end{array}\right) . \end{aligned}$$

(5)

According to Refs. [35, 49], this Hamiltonian matrix includes the effective Hamiltonians ${H}_{\mathbf{k}{\mathbf{k}}^{\mathbf{\text{'}}}}^{\alpha \alpha \text{'}}\left(\mathbf{K}\right)$ , where $\mathbf{k},{\mathbf{k}}'$ are the wave vector of an electron in the Brillouin zone of the crystal, $\alpha ,\alpha'$ are the indexes of the considered bands (${\Gamma }_{6}^{c}$, ${\Gamma }_{8}^{v}$ and ${\Gamma }_{7}^{v}$ are the centers of the Brillouin zone of a bulk crystal for the conduction, valence, and spin-split bands, respectively, with spin momentum projections ${\Gamma }_{6}^{c}(c= \pm 1/2)$, ${\Gamma }_{8}^{v}(v=\pm 3/2, \pm 1/2)$ and ${\Gamma }_{7}^{v}(v= \pm 1/2)$, $\mathbf{K}$ is the wave vector of an electron in the Brillouin zone of SL. The main diagonal elements are the one-band effective Hamiltonians, and the off-diagonal elements describe the effects of the interband kp-interaction. Both types of Hamiltonians include the difference in the interband parameters of the materials when passing through the interface. When solving this problem, one can prevent false solutions in the multiband model by the $\mathbf{k}$—representation method since the considered region of $\mathbf{k}$—space can be limited based on the estimated width of the wave functions compared to the size of the Brillouin zone.

The primary method to solve the Schrödinger equation involved a periodical continuation of a single quantum well^{[35, 51, 52]}. The periodicity condition leads to the fact that the wave vector $\mathbf{k}$ of an electron in a crystal becomes discrete, and the Schrödinger equation itself turns into a system of algebraic equations

$$\sum _{{\alpha }',{\mathbf{k}}'}{H}_{\mathbf{k}{\mathbf{k}}^{\mathbf{\text{'}}}}^{\alpha {\alpha }'}\left(\mathbf{K}\right){F}_{{\mathbf{k}}^{\mathbf{\text{'}}},l}^{{\alpha }'}\left(\mathbf{K}\right)={E}_{l}\left(\mathbf{K}\right){F}_{\mathbf{k},l}^{\alpha }\left(\mathbf{K}\right),$$

(6)

where ${F}_{\mathbf{k},l}^{\alpha }\left(\mathbf{K}\right)$—the envelope wave function of an electron in the $\bf{k}$-representation, $l$ is the number of the band in the SL band structure. The components of the eigencolumns of the Hamiltonian ${F}_{\mathbf{k},l}^{\alpha }\left(\mathbf{K}\right)$ determine the contribution of the states $\left|\alpha ,\mathbf{k}\rangle\right.$ of the reference crystal corresponding to the bands to the superlattice state $\left|l,\mathbf{k}\rangle\right.$.

To interpret the experimental data presented in Fig. 1, the calculation of the energy spectrum of charge carriers by Eq. (6) is not enough, since the experimental data reflect only the optical transition probabilities. The optical transition probabilities can be obtained from analyses of interband matrix elements, which are calculated with the following formula (neglecting the kp-interaction of electrons and holes):

$${\mathbf{P}}_{l{l}'}\left(\mathbf{K}\right)=\sum _{c,v}\sum _{\mathbf{k}}{{F}_{l,\mathbf{k}}^{{\Gamma }_{6}^{c}}}^{\dagger}\left(\mathbf{K}\right)\cdot {\mathbf{P}}_{cv}{\cdot F}_{\mathbf{k},{l}'}^{{\Gamma }_{\mathrm{8,7}}^{v}}\left(\mathbf{K}\right),$$

(7)

where ${\mathbf{P}}_{cv}$ is the matrix element of the momentum operator. For transitions of charge carriers from the valence band ${\Gamma }_{8}^{v}$ to the conduction band ${\Gamma }_{6}^{c}$, ${\mathbf{P}}_{cv}=\langle{\Gamma }_{6}^{c}\left|\widehat{\mathbf{P}}\right|{\Gamma }_{8}^{v}\rangle$. For transitions from the spin-split band ${\Gamma }_{7}^{v}$ to the conduction band ${\Gamma }_{6}^{c}$ ${\mathbf{P}}_{c{v}_{so}}=\langle{\Gamma }_{6}^{c}\left|\widehat{\mathbf{P}}\right|{\Gamma }_{7}^{v}\rangle$, its projections have the form

$$\begin{split} {P}_{cv}^{x}=P\left(\begin{array}{cc}\dfrac{-i}{\sqrt{6}}& 0\\ 0& \dfrac{-i}{\sqrt{3}\sqrt{6}}\end{array}\begin{array}{cc}\dfrac{i}{\sqrt{3}\sqrt{6}}& 0\\ 0& \dfrac{i}{\sqrt{6}}\end{array}\right),\\{P}_{cv}^{y}=P\left(\begin{array}{cc}\dfrac{1}{\sqrt{6}}& 0\\ 0& \dfrac{1}{\sqrt{3}\sqrt{6}}\end{array}\begin{array}{cc}\dfrac{1}{\sqrt{3}\sqrt{6}}& 0\\ 0& \dfrac{1}{\sqrt{6}}\end{array}\right),\\{P}_{cv}^{z}=P\left(\begin{array}{cc}0& \dfrac{i}{\sqrt{3}\sqrt{6}}\\ 0& 0\end{array}\begin{array}{cc}0& 0\\ \dfrac{i}{\sqrt{3}\sqrt{6}}& 0\end{array}\right),\quad\;\; \end{split}$$

$$\begin{split}& {P}_{{cv}_{so}}^{x}=P\left(\begin{array}{cc}0& \dfrac{i}{3}\\ \dfrac{i}{3}& 0\end{array}\right),\\&{P}_{{cv}_{so}}^{y}=P\left(\begin{array}{cc}0& \dfrac{i}{3}\\ \dfrac{-i}{3}& 0\end{array}\right),\\&{P}_{{cv}_{so}}^{z}=P\left(\begin{array}{cc}\dfrac{i}{3}& 0\\ 0& \dfrac{-i}{3}\end{array}\right),\end{split}$$

(8)

where $P=\langle{\Gamma }_{1}\left|{\widehat{p}}_{x}\right|{\Gamma }_{15}^{x}\rangle=\langle{\Gamma }_{1}\left|{\widehat{p}}_{y}\right|{\Gamma }_{15}^{y}\rangle=\langle{\Gamma }_{1}\left|{\widehat{p}}_{z}\right|{\Gamma }_{15}^{z}\rangle.$ Since the structures under study were grown in the direction [001], the projection of the matrix element onto the z axis did not contribute. It is important to note that the band anticrossing effects take place in the band structure for the hole states of the superlattice in the ${K}_{z}$ direction. This leads to a sharp change in the selection rules for the optical matrix element ${\mathbf{P}}_{l{l}'}\left(\mathbf{K}\right)$ when passing through the ${K}_{z}$ point corresponding to the "anticrossing" of terms.

For the energy spectrum of charge carriers and interband matrix elements calculations the parameters of initial materials from Table 2 were used^[54]. The Kane parameters required for the calculation were determined through the constants presented in Table 2 as follows:

$${\gamma }_{c}^{K}={\gamma }_{c}^{L}-\frac{{E}_\rm{p}}{3}\left[\frac{2}{{E}_\rm{g}^{\left({\Gamma}\right)}}-\frac{1}{{E}_\rm{g}^{\left({\Gamma}\right)}+{\Delta }_\rm{so}}\right],{\gamma }_{c}^{L}=\frac{1}{{m}_\rm{e}^{*}\left({\Gamma}\right)},$$

$${\gamma }_{1}^{K}={\gamma }_{1}^{L}-\frac{{E}_\rm{p}}{3{E}_\rm{g}^{\left({\Gamma}\right)}},{\gamma }_{2}^{K}={\gamma }_{2}^{L}-\frac{{E}_\rm{p}}{6{E}_\rm{g}^{\left({\Gamma}\right)}},{\gamma }_{3}^{K}={\gamma }_{3}^{L}-\frac{{E}_\rm{p}}{6{E}_\rm{g}^{\left({\Gamma}\right)}},$$

$${\gamma_{i}^K}_{\mathrm{Al}_x\mathrm{Ga}_{1-x}\mathrm{As}}={\gamma_{i}^K}_{\mathrm{AlAs}}\cdot x+{\gamma_{i}^K}_{\mathrm{GaAs}}\cdot\left(1-x\right),$$

$${E}_{p}=\frac{2{P}^{2}}{{m}_{0}},$$

(9)

where ${m}_{0}$ is the electron mass.

The band offset parameters were obtained from Ref. [50] and corresponded to the following relationship ${\mathrm{\Delta }E}_\mathrm{c}/\mathrm{\Delta }{E}_\mathrm{v}= 0.67/0.33$. The band offset of spin-split bands was defined as $\mathrm{\Delta }{E}_\rm{SO}=\mathrm{\Delta }{E}_\mathrm{v}+\left({\mathrm{\Delta }}_\rm{SO}^\mathrm{GaAs}-{\mathrm{\Delta }}_\rm{SO}^{\mathrm{A{l}}_{x}\mathrm{G{a}}_{1-x}\mathrm{As}}\right)$. It was also assumed that neglected deformation effects are insignificant since the relative difference between the constant lattices of the well and barrier materials is $\left|{\epsilon }_{\perp }\right|=0.097 \%$. Therefore, the lattice period in the superlattice was taken as the Al_xGa_1-xAs period. We computed the transition energies for considered SMP SL structures using the above-described formulae. The resulting values are demonstrated in Table 3 and they were used for further PR spectra analysis.

3.3 Theoretical calculations results vs experimental spectra

The above-mentioned results can be also illustrated as the electronic minibands in the conduction band in Fig. 2. Figs. 2(a) and 2(b) show the relevant E(K) dependences for the electronic states, where K varies within the Brillouin zone of the SL ($(K_{_z}\neq0,K_{_x},K_{_y}=0)$). The sizes along the z-dimension of the unit cell of the SLs in the model are D = 12.4 and D = 9.05 nm, respectively. From these figures, one can conclude that the third electronic subband e3 for SL no. 1 is partially localized within the well. As a result, its PR spectrum has optical transitions involving e3 not only at point Г but also at point П (edge of Brillouin zone of the SL) with above-barrier energy levels (see Fig. 2(с)).

Within the scope of the considered, it was assumed that the profile of quantum wells in the superlattice is rectangular (trapezoidal—under electric fields) with abrupt interfaces. The distortions of the conduction/valence bands profiles were neglected due to the insignificant magnitude for the considered doping concentrations. The presence of sharp interfaces was shown by the results of HR HAADF STEM, HRXRR, DXRR, HRXRD investigations. To compensate the minor distortions of the potential profile, the size of a single SL period for the first sample (with period 12 nm) was increased by one GaAs monolayer to obtain better convergence.

Electronic energy spectra calculated by the Kane model (Table 3, Fig. 2) were consistent with transformed PR spectra of the synthesized Al_0.3Ga_0.7As/GaAs SMP SLs (Fig. 3). The comparison between experimental and theoretical data allowed us to interpret transformed PR spectra in detail. All the simulated transitions are presented on obtained spectra in Fig. 3 with colored solid arrows. The small broadening of spectral lines shown in Fig. 3(a) confirms the high quality of SMP SL samples. The most intensive and narrow spectral line at 1.449 eV corresponds to the main optical interband transition e1−hh1 at the Г point (center of the Brillouin zone of the SL). The next less intensive peak at 1.461 eV corresponds to the e1−lh1 transition at the Г point. The high quality of the SL has led to a small broadening of spectral lines and these two peaks with extrema position difference of 12 meV turned out to be resolved even at room temperature. Their position, as well as the position of other marked peaks, is in good agreement with the Kane model calculations.

Fig. 3 also compares the transformed PR spectra of the superlattices samples no. 1 (Fig. 3(a)) and no. 2 (Fig. 3(b)) with periods D = 12.0 nm and 8.99 nm, respectively. A decrease in the superlattice period from 12.0 to 9 nm led to a significant increase in the same optical transition energies observed in the spectrum (Fig. 3(b)). The main optical transition in sample no. 2 increased in energy up to 1.483 eV. The left shoulder of the main optical transition peak in the transformed spectrum is assumed to be associated with the impurity exciton state. According to Kane model calculations with such SL configuration, two electron subbands are formed, while the second subband is only partially localized within the well (see Fig. 2(с)). As a result, relatively intense transitions are observed in the spectrum with the participation of the e2 subband at the П point around 1.63−1.68 eV and less intense transitions at the Г point with above-barrier transitions.

The difference between experimentally obtained and computed transition energies does not exceed several meV. The proposed Kane model allowed us to significantly reduce the error of the peak positions from 25% to 1%−5%. Also, using the Kane model, it became possible to obtain the precise position of minibands in the range of SL designs (Fig. 4). Obtained results show smooth dependence from the width of the well. It is also worth mentioning that there is a threshold width value of 5 nm. For more narrow quantum wells, the expectation is that the second eigenstate won’t be observable in PR analysis.

The Kane model described above is highly accurate, but it can simulate the band diagram and electronic energy spectrum only. This model is not intended to calculate electronic transition rates between corresponding energy states. Therefore, other methods, that consider electronic transitions, should be used to estimate scattering rates and gain. Since CSM, described earlier, can calculate these characteristics, we decided to combine the 8-band kp-model with the “shooting approach”. The simulation process of our method can be divided into the following steps:

(1) initially, energy level positions (their vicinity), effective masses, theoretical period of SL, etc. were computed with the Kane model;

(2) computed parameters were fed to CSM to obtain the accurate solution of the Schrödinger equation, i.e., the eigenvalues and eigenvectors for the given Hamiltonian;

(3) the figures of merit, transition rates, and gain estimates were calculated by MC approach and Eqs. (1)−(4) with eigenvalues and eigenvectors, obtained in the previous step.

This sequence of steps preserves the accuracy of the Kane model and allows to estimate gain. Thus, this method has all the prerequisites to become a new and more reliable alternative to commonly used approaches.

To get the highest gain estimate and maintain the frequency of amplified light in the THz range, we have varied the thickness of the quantum well between 7 and 10 nm. The results obtained by our new method are presented in Fig. 5−Fig. 7. All of the below-mentioned computations were performed for the temperature of 300 K. Fig. 5 demonstrates the frequency dependence of the estimated gain for different quantum well thicknesses. The graph illustrates that wider quantum wells provide higher gain values, but the radiation frequency is blue-shifted towards 2.5 THz. Presumably, the higher gain values are achieved due to the increased wave function localization in wider wells.

Fig. 6 shows the figure of merit dependence on the quantum well thickness. The plot confirms that wider quantum wells result in higher gain, as shown by the figure of merit dependence on quantum well width. According to this graph, the most suitable width values, where the highest F_M can be obtained, are between 8.5 and 9.5 nm. The maximum position is expected to be at 8.8 nm.

Scattering rates for LO-phonon and impurity mechanisms were also calculated during simulations, and their dependency on quantum well width is shown in Fig. 7. Shorter scattering times result in faster depopulation of the level, i.e., designs with longer scattering times for the upper level and shorter times for the lower level are the most convenient for the population inversion maintenance. Therefore, the most suitable region in Fig. 7 is between 8.5 and 10 nm:

• scattering times for the lower level at their minima (they do not exceed 1.5 ps);

• scattering times for the upper level are relatively high (approximately 5 ps).

To summarize all the mentioned results, the most promising structure for a THz source implementation is a considered SMP SL with an 8.5 nm thick quantum well. This geometry provides a higher gain estimate and allows to maintain population inversion with small broadening values $\gamma$.

The described combination of PR and Kane model can be used in a more self-contained manner. To show that, the authors also enhanced this method with the process of naïve optimization, where all the following ranges of SL parameters are brute-forced: mole fraction of Al in barrier layers between 10% and 40%, well thickness from 4 to 15 nm. Achieved minimum deviation met expectations, and its position was in good agreement with XRD data and growth parameters. By changing the thickness of the well, one can also control the number of eigenstates inside the SL. As was previously shown, wider quantum wells contain a third undesirable energy level^[9]. The proposed method predicted the escape of the second energy level when reducing quantum well thickness, which is illustrated in Fig. 2.

To minimize the adverse impact of elastic scattering, it is suggested to increase the number of periods and decrease the doping concentration of the layers. According to the estimations (Fig. 5), this must significantly decrease (down to several magnitudes) the influence of the elastic scattering, while maintaining high gain values. With the above-described adjusted energy values from the Kane model, the dipole matrix element, relaxation times, and the estimation of the gain coefficient were calculated (Fig. 5−Fig. 7). Obtained values showed a gain slightly above 40 cm⁻¹ for the proposed SL with a quantum well layer thickness of 8 nm. We assume that the gain values (see Fig. 5) are significant enough to produce light emission, as the total expected losses are estimated to be below 15 cm⁻¹. To estimate the waveguide loss, we considered the double metal waveguide from Ref. [55]. It should be also mentioned that this structure is the best for light amplification due to the electron transport estimations that exhibit a positive differential conductivity region near the operating point^{[9, 11]}.

4. Conclusion

In conclusion, this paper presents the first and unique comprehensive theoretical analysis of the energy band diagram for super-multiperiod superlattices, verified by experimental data. The authors have developed a new combined method for the theoretical study of multiperiod superlattices. It fuses the precision of the Kane model and the flexibility of the shooting method with the Monte Carlo technique, which results in a more refined and versatile approach. This method was applied to analyze state-of-the-art samples of super-multiperiod superlattices grown by molecular beam epitaxy. The express room-temperature method of photoreflectance was used for experimental approbation. The proposed method diminished the error in the identification of the transition energy peaks to just 5%, presenting a substantial improvement compared to the typical 25% error obtained with conventional approaches. The photoreflectance experiments validated the accuracy and reliability of the combined method. The approach successfully predicted the escape of the second energy level in structures with quantum well thicknesses below 5 nm. The authors reinforced the method with naïve optimization to restore the intrinsic superlattice parameters. The recovered parameters correlated well with X-ray diffractometric data and growth parameters.

The refined method is also suitable for calculating transition rates and estimating gain. Therefore, it was used to design a new light-amplifying super-multiperiod superlattice. Additionally, the authors revealed ways to reduce the impact of elastic scattering in these constructions. The scheme developed in the research boosts the creation of both efficient light amplification and adjustable sources of terahertz and infrared radiation.

Appendixes A and B. Supplementary materials

Supplementary materials to this article can be found online at https://doi.org/10.1088/1674-4926/45/2/022701.https://doi.org/10.1088/1674-4926/45/2/022701.

Acknowledgements

The work was supported by the Ministry of Education and Science of the Russian Federation in the framework of experimental research (Nos. 075-01438-22-06 and FSEE-2022-0018) and the Russian Science Foundation in theoretical research (No. RSF 23-29-00216).