Dual resonance phonon–photon–phonon terahertz quantum-cascade laser: physics of the electron transport and temperature performance optimization

: The state of the art terahertz-frequency quantum cascade lasers have opened a plethora of applications over the past two decades by testing several designs up to the very limit of operating temperature, optical power and lasing frequency performance. The temperature degradation mechanisms have long been under the debate for limiting the operation up to 210 K in pulsed operation in the GaAs/AlGaAs material system. In this work, we review the existing designs and exploit two main temperature degradation mechanisms by presenting a design in which they both prove beneficial to the lasing operation by dual pumping and dual extracting lasing levels. We have applied the density matrix transport model to select potential candidate structures by simulating over two million active region designs. We present several designs which offer better performance than the current record structure.


Introduction
Quantum cascade lasers (QCLs) have been initially proposed by Kazarinov and Suris in 1971 [1] as periodic two level resonant tunnelling superlattice structures. The first realization came in 1994 [2] at 4.2 µm which gave rise to multiple designs that lase across mid-infrared spectrum with high power and room temperature operation [3,4]. The paradigm of maintaining population inversion via resonant tunnelling has shown promise for scaling towards higher frequencies and reaching the traditionally "invisible" far-infrared portion of electromagnetic spectrum. In 2002 the first terahertz-frequency (THz) QCL [5] has been realized, lasing at 4.4 THz up to 50 K in pulsed operation. Since then, a variety of designs [6] generated structures with lasing frequencies in the range 1.2-5.6 THz [7][8][9][10], with high output power [11] and with operating temperature up to 210 K [12,13] in pulsed and 129 K [14] in continuous wave (CW) operation (without the assistance of external magnetic field). These limits have all been achieved in GaAs/AlGaAs heterojunction superlattice which has shown the best performance for THz technology, however the hope for improvement currently lies in other material systems [15,16] which exhibit better material parameters, but require further technology improvements.
The lower frequency limit of 1.2 THz is expectedly present since the required lasing levels energy separation is just ∼4.9 meV, while the upper frequency limit is a consequence of highly absorbent Reststrahlen band [17].
The temperature limit comes from the very need of small lasing levels energy separation (of ∼12 meV), detrimental scattering of electrons and longitudinal optical (LO) phonons which in GaAs is most dominant at ∼36 meV, and parasitic leakage between states which do not participate in lasing process directly, but cannot be avoided by the design itself. All these problems are hard to avoid and solutions to particular issues usually counteract the other ones, leading to compromising designs. The best THz QCL designs for high temperature performance [12,13] are usually focused around 3.2-3.8 THz where the waveguide loss is predicted to be minimal [18], the LO-phonon scattering is exploited to assist the lasing process, and parasitic leakage may be suppressed by using higher barriers [13] which on the other hand causes higher thermal heating of the device.
Historically, vast majority of lasers may be represented by three-level and four-level schemes [19]. In Section 2, we classify common THz QCL designs by effective three-and four-level schematics, and discuss their performance issues at high temperature. In Section 3 we formulate a novel four-level design which takes advantage of detrimental processes to act beneficially towards maintaining the population inversion. Section 4 proposes several structures that match our design paradigm. We used our numerically highly efficient density matrix transport model [20,21] that allowed brute-force search by simulating over two million candidate structures.

THz QCL designs
Lasing process in THz QCL is based on the same principles as in almost all other semiconductor, solid or gas lasers: lasing transition occurs between two levels -upper lasing level (ULL) and lower lasing level (LLL) where the main aim is to maintain population inversion by efficiently injecting carriers into ULL and extracting them from LLL. In many cases intermediate lasing level (ILL) is introduced to assist the pumping/extraction process. Every QCL exploits resonant tunnelling effect usually to extract LLL and pump ULL or ILL in the next period. Multiple ILLs may be introduced to form three, four or more level lasing systems and it is also possible to form LLL and/or ILL as a miniband with narrowly spaced states to achieve faster relaxation of the levels.
Band-structure designs of THz QCLs are usually classified as bound-to-continuum (BTC), resonant (or LO)-phonon, and hybrid designs [22]. The latter combines advantages of the former two, however resonant phonon structures display the best temperature robustness [12,13]. An additional classification introduces scattering-assisted designs [23][24][25][26] in which the role of pumping/extraction of lasing levels processes is either reversed (in comparison to typical LO-phonon design) or enhanced (ILLs are added to the design so that both pumping and extraction are performed via LO-phonon scattering).
If we adopt a simplified view where we allow any state apart from ULL to be formed out of a mini-band of narrowly spaced levels, we can represent the majority of THz QCL designs as effectively three-and four-level lasing systems, depicted in Fig. 1 [13,[23][24][25][26][27][28][29].
The simplified schematics shown in Fig. 1 may be applied to nearly every experimentally realised THz QCL. The hybrid design would simply represent the case of resonant phonon structure in this schematics, where LLL and/or ILL are minibands instead of discrete levels. Note that original classification to BTC, LO-phonon, scattering assisted and hybrid structures is not definite and in the literature a hybrid design would often be referred to as scattering-assisted BTC design, or scattering-assisted designs would all refer to schemes in Figs. 1(c)-(f).
QCL is a periodic structure, and resonant tunnelling occurs at the resonant electric field bias which causes alignment of states from two adjacent periods in such a way that one period performs the pumping process to the second period, while the second period performs the extraction process from the first period or vice versa. The resonant bias K across QCL period of length L P therefore depends on the required energy splitting ∆E R needed to bring the desired states into Fig. 1. Schematic diagram of effectively two (a), three (b,c) and four (d,e,f) level schemes of common THz QCL designs. The rectangles illustrate the typical wavefunction localisation (probability density) of each state within the QCL period. The dotted arrow line denotes the tunnelling process between two adjacent periods, while the solid arrow lines denote the transitions. Each level (apart from ULL) may be envisaged as a cluster of narrowly spaced quasi-bound levels, transitions between the effective "levels" also exist, however dominant mechanisms are shown. resonant tunnelling position: where ω is the lasing frequency, ℏω LO is LO-phonon energy (∼36 meV in GaAs), ∆E minibands is added to account for cases in which some of the states in effective schematics in Fig. 1 consist of a cluster of narrowly spaced states, and indices a-f refer to particular values of ∆E R for cases presented in Fig. 1.
The BTC design in schematics in Fig. 1 uses tunnelling effect to pump ULL in the next period, however the extraction from LLL within the period needs a miniband which does not lead to good performance at higher temperatures. The electron-LO-phonon scattering mechanism is highly dominant in semiconductor materials and although originally detrimental, all designs (apart from BTC) exploit its efficiency to implement some form of ILL. The following processes have detrimental effect on the temperature performance [26]: 1. Electron-LO-phonon scattering activates on a similar energy scale as the lasing energy difference at higher temperatures. It has a peak around 36 meV in GaAs, however the effect persists for a wide range of energy differences (∼ 12 -60 meV at high temperature [30]). If this mechanism is used to assist extraction of LLL (typical designs place ILL 36 meV below LLL), it may also inherently extract carriers from the ULL (i.e. in typical resonant phonon design ULL may be ∼15 meV above LLL). Similarly, if this mechanism is used to pump ULL (i.e. typical scattering-assisted design) it may also pump the LLL. We will refer to this detrimental effect as LO-phonon leakage.
2. Tunnelling effect requires the two states from adjacent periods to about hybridize and reach resonant energy difference equal to their anticrossing energy. It is clear that tunnelling process may also interact with an undesired level. We will refer to this effect as tunnelling leakage.
3. Electrical heating dictates many trade-offs in THz QCL designs. The electrical power depends on current and applied voltage (which is proportional to the bias). Resonance bias depends on design, as illustrated in Eq. (1), and is dependent on the period length and the desired energy difference. Since LO-phonon energy in GaAs is ∼36 meV, the designs in Fig. 1 that exploit two LO-phonon scattering processes need higher bias. 4. Period length has a twofold effect on temperature performance: it affects electrical heating and it also affects the number of states per period. This causes a trade off, because if a miniband of states is formed in THz QCL, the temperature performance drops down, since a narrow cluster of states may exhibit multiple undesired absorption and non-radiative processes at high temperature. On the other hand, short period designs have more pronounced dependence on ∆E R , which leads to higher values of current and voltage threshold, and thus increases the electrical heating. For this reason, BTC and hybrid THz QCL designs perform well in CW operation, while typical resonant phonon devices operate poorly in CW, but are highly suited for high temperature performance in pulsed operation.
5. Parasitic leakage may occur due to quasi-bound levels that are not part of the desired transport schematics (or are in continuum), but cannot be avoided. For instance, construction of a three level resonant phonon THz QCL can be obtained by using only two quantum wells [13]. One well needs to be very wide to generate ILL and LLL, while the second one needs to be narrow in order to form ULL at a high energy. Unfortunately, the wide well will also cause formation of additional higher levels that may trigger thermal backfilling. This can be suppressed by using higher barriers [13], however this would increase the effective mass (thus reducing gain), and require higher doping which will cause higher electrical heating.
Analyzing the temperature performance of all designs in Fig. 1 would reveal a different trade-off between the processes mentioned above. The current record temperature performance limit is 210 K [13]. This structure implements a two-well design that directly generates the three quasi-bound states. However, two-well designs need high barriers in order to suppress the parasitic leakage, and this design has very high doping, and therefore very high threshold (24 V, 9.5 A), which causes significant electrical heating, that was avoided by using very low duty cycle in [13]. We can formulate a rough electrical heating prediction [30] which uses a linear dependence T = T 0 + δR TH IV, where T 0 is the heat sink temperature, R TH is the thermal constant, δ is the duty cycle, I is the current and V the voltage. Note that transport models usually deal with current density J and electric field bias K, which can be approximately taken to scale uniformly across the emitting surface and active region ridge height, respectively. We can describe the electrical heating by introducing the heating factor, as follows: ]︂ , meaning that a QCL driven with 2% duty cycle would have ∼ 32 K higher temperature than the cold finger temperature. Therefore, any prediction of material gain dependence on temperature would show a significant shift. Note that in [13] the duty cycle was 0.01 % in order to avoid this effect.
We will also focus specifically on transport leakage effects in designs presented in Fig. 1. The tunnelling leakage was often discussed as the main limiting process in the temperature performance. Several "exotic" designs with variable barrier heights [31] have been attempted, and scattering-assisted structures have shown that tunnelling effect may be used for LLL extraction instead of ULL pumping. Although several theoretical works [32][33][34] propose three-level scattering-assisted THz QCLs with shorter injection barriers, nearly all realised devices still use the barrier thicknesses as in the resonant-phonon devices. This may be attributed to tunnelling leakage, or to a lack of experimental design exploration, however the main issue of these structures is typically a high operating bias (and therefore electrical heating). A design that does not experience any tunnelling leakage to undesired levels is phonon-photon-phonon (PPP) structure [23,24,26], and a high operating temperature was achieved at low frequency [23,26], however the design for a higher frequency [24] did not deliver any significant improvement despite lower waveguide loss. Performance of PPP design may be attributed to the need for a significantly higher operating bias and the leakage caused by LO-phonon scattering process.
It is clear that LO-phonon scattering leakage mechanism may appear in all the presented threeand four-level schemes. In PPP design this may be twice worse than in the resonant-phonon structure since both ILL states may start pumping/extracting the undesired lasing level. Prevention of LO-phonon leakage may be achieved by designing ILL states to be >ℏω LO away from ULL/LLL. This is the case in [13], and will be the aim of the design that we will present in the next section.

Dual resonance phonon-photon-phonon THz QCL design
The following design has been inspired by the three-level scattering-assisted design [25] and PPP structure [23,24,26]. It is rather clear that the issue with the three-level scheme in [25] was in inefficiency of tunnelling to extract carriers from the LLL, and while the PPP design has solved this issue, it caused additional problems by operating at very high voltage and also led to increased LO-phonon leakage. Note that a large number of PPP structures have been realized in [6], where some levels were formed as minibands, however this did not provide any better temperature performance, as expected. Interestingly, a variant of the PPP design presented in Fig. 2 has (to our knowledge) never been attempted or proposed in the literature so far. The structure proposed in Fig. 2 can be viewed as a four-level phonon-photon-phonon design that uses two tunnelling processes, or as a three-level resonant phonon scheme where ILL 1 is added to act as a parasitic state that assists the extraction and pumping processes. In either case, this proposed design offers the following advantages: • The required potential drop is identical to that in the three-level resonant phonon designs (in Eq. (1): ∆E R = ℏω + ℏω LO .
• Pumping mechanism of ULL is achieved by the usual tunnelling process as in the resonant phonon design and the additional LO phonon process with a state above which is in resonance with LLL from the previous period.
• Extraction of carriers from LLL is achieved by the usual LO phonon process as in the resonant phonon design and the additional resonant tunnelling process with a state from the next period, that would undergo efficient LO-phonon scattering process with ULL in the adjacent period.
The potential disadvantages are: • The design may require high barriers to generate the ILL 1 state in order to suppress the higher parasitic states created by the widest well (that is needed to create ILL 2 and LLL) and to achieve operating frequency that has low loss (3.2 -4 THz). Higher barriers may also be needed to suppress the leakage into continuum, since ILL 1 is a high state.
• The dual tunnelling process may be affected by growth, as this design needs precise layer widths to avoid parasitic depopulation of ULL.
This design has high symmetry in its transport mechanisms. It essentially implements threelevel scattering assisted design over resonant phonon design, with the aim to collect and recycle any potential leakage. This design seemingly solves the parasitic processes issues that may be responsible for temperature degradation in designs discussed in the previous section. The full symmetry of transitions acts in the following way at high temperature: if ILL 1 starts to pump the LLL as well, a part of population of LLL will be re-injected into ILL 1 of the next period, and similarly if ILL 2 starts to extract the ULL as well, a part of population of ULL will be injected into ULL in the next period. Similar issues arise in case that one of the tunnelling processes starts to pump/extract the undesired state.
This design paradigm has been obtained through the review presented in the previous section. The next challenge is to find the layer sequence which matches the design criteria. Note that potential parasitic states occurring in the vicinity of ILL 1 would be highly detrimental. Additionally, designs similar to this one do not exist in the literature (to the best of our knowledge), thus it is not possible to just tune the layers' widths around the values in already known designs.

Numerical results
In the previous work [20,21] we have presented a density matrix model applicable for any number of states per module, and found that it is capable of very good estimation of the cut-off temperature for a variety of QCL designs [35]. This DM model is a generalization of the model from [36] that was used to optimize the earlier record high operation temperature structure [12]. Another model, very similar to ours, has also shown good agreement [37] in the cut-off temperature prediction. For instance, the theoretical model for 210 K record structure [38] predicted positive material gain over 20 cm −1 up to 300 K, and waveguide loss around 30-35 cm −1 , while our density matrix model [20,21] predicts the material gain of >30 cm −1 exactly up to 210 K and pronounced deterioration of material gain towards higher temperatures. Similarly, our model predicts exactly the cut off temperature of 200 K [12] structure for which the authors predicted the loss around 17 cm −1 .
Our DM model, [20,21] uses tight binding approximation for determination of quasi-bound states in one QCL period meaning that we set the QCL period layer structure starting and ending with half a thickness of the injection barrier. Then the injection barrier halves are additionally extended to infinity on both ends (numerically this is satisfied by using padding around the period ends of ∼ 20 nm). We define the injection barrier to be the barrier between injector state ILL 2 and upper lasing state ULL from neighbouring period in Fig. 2. The Schrödinger -Poisson equation is then solved in self-self-consistent manner where we use equithermal approximation for subband temperatures which determines the electron temperature through minimization algorithm. Such solutions are then plugged into DM model that requires information on scattering mechanisms within periods and coherent effects between period's nearest neighbours. The scattering effects we considered in this work include interactions of electrons with optical and acoustic phonons, alloy disorder, interface roughness and ionised impurities. Carrier-carrier scattering has not been included due to its heavy computational load. The inter-period coherent interaction in our DM model requires information on Rabi coupling strengths energies (proportional to anticrossing energy if states are in resonance) and we used a model in [39]. A more general approach is available in [40]. Note that since we are extending the injection barrier due to tight binding nature of the model, this barrier does not affect the energies and wavefunctions in our model, however the actual injection barrier thickness affects the calculation of Rabi coupling strengths.
Non-equilibrium-Green-Function (NEGF) models are frequently used in the literature [13,33], however they typically suffer from very high computational cost. The algebraic simplification we have made in [21] allowed simulation of three-and four-well structures in 8-15 s per single bias input, where the simulation time bottleneck is not the DM model, but rather the Schrödinger-Poisson solver. The typical approach for QCL optimization would be to fix the operating bias in accordance to the desired transition, as in Eq. (1), and vary the critical layers of the structure [38]. However, with the high speed performance of our model and access to high performance computing cluster (ARC4), we were able to vary all layers of three-and four-well QCL structures and sweep the bias over 40 steps while keeping the barrier heights, doping and injection barrier widths constant, thus performing a brute-force numerical approach.
Simulation of several million THz QCL structures was performed in searching for the best designs that meet the paradigm presented in Fig. 2. The initial methodology was focused on three-well THz QCL, where the aim was to take a two-well THz QCL and add one narrow well to generate ILL 1 state at high energy (>36 meV above ULL). The findings of these simulations pointed that a four-level design would perform better, and we then conducted an even more extensive set of simulations for four-well structures. We should note that all layers apart from the injection barrier were varied, since our model is tight-binding, and the lattice temperature was fixed at 250 K.
Although our design shares the properties of both the resonant phonon and scattering-assisted design, we decided to follow the common doping practice by doping only the widest well, as in the resonant phonon design. We chose to dope mostly the second half of the well. The common procedure in latest two temperature records [12,13] was to dope the narrow central portion of the widest well, however in the present design we assumed that it is not detrimental to shift doping towards the injection barrier, as we want another resonant tunnelling process. It should be noted that the output of our model shows negligible dependence on what part of the well is doped, however this would need to be taken into account in the experimental growth.
There is some inconsistency between theory and experiment in literature regarding the injection barrier width in scattering-assisted designs. As discussed in Section 2, there are several proposals with thinner injection barriers [32][33][34], while the realised structures [23,24,26] (that we found) used the injection barriers thickness typical for LO-phonon structures. For this reason we kept the values commonly used in LO-phonon designs, however from most simulations it appears that using a thinner injection barrier may yield higher material gain, as will be discussed.

Three-well design
The layer variation details for three-well design are presented in Table 1 where we simulated ∼500,000 candidates.
The simulations in Table 1 were performed for Al molar fractions x = 0.15, 0.16, 0.17, 0.23, 0.24, 0.25, where we changed the injection barrier thickness to 42.375, 42.375, 39.55, 33.9, 33.9, 33.9 Å, respectively, as is commonly done in designs with high barriers. Only the widest well has been doped. 84.75 Å of the well was undoped, while the rest of the well's width was doped so to correspond to 2 · 10 10 cm −2 sheet doping density for all Al fractions. The lattice temperature was 250 K.
The focus on barriers with x = 0.15 − 0.17 is because x = 0.15 is a value typical in majority of THz QCLs, and it is also desirable in most experimental equipment calibration. Additionally, the higher barriers in QCL design have higher effective mass and therefore may have lower material gain for the same value of sheet doping density. The structures with high value of the material gain did in fact follow the paradigm presented in Fig. 2, however we noticed that state ILL 1 (the highest state in the period) would also be the highest state set in the simulation. In a three well design, we are interested in finding layer sequence with four levels in a period and hope that higher states would not be too detrimental to transport. However, if e.g. six states are allowed in the simulation, the best candidates would have state six as ILL 1 and we noticed that if we add just one more state above, the material gain deteriorates because state seven would most likely be in close resonance with ULL or cause some other form of leakage. The initial paradigm in Fig. 2 already requires ILL 1 to be quite high, and any state higher than that would be nearer to the continuum, thus we have decided that this set of simulations that predicted high material gain with six and seven states (the eight state is typically in the continuum) are invalid.
This issue was rectified by implementing the following algorithm: • Perform simulation with only four states included • Detect the subset of simulations with high value of material gain (> 15 cm −1 ) • Perform additional simulations for the selected subset by allowing five states, six states and seven states per period • The final analysis is focused on the data obtained from the seven-state simulation, with caution, since state seven may have ended up in the continuum.
In this way we were investigating the effect of leakage caused by undesired higher states present in the simulation. Ideally, structures that do not suffer from significant material gain or lasing frequency change with the number of states in the simulation would be the best candidates, and we have indeed found such structures only for higher barriers.
Interestingly, simulations done for x = 0.18 − 0.22 (not presented in Table 1) generated candidate structures with low material gain and with slightly different transport principle than the one proposed in Fig. 2: The state above ILL 1 would be in resonance with ULL from the previous period, thus creating a dual lasing channel. However this channel would not be very efficient because this high state has no efficient way of being pumped and acts mainly as a parasitic state.
We introduce a nomenclature for our designs as follows: D_3_x_# where 3 means that this is a three-well structure, x (in [%]) is the Al mole fraction and # is the simulation number obtained through nested loop variation of the ranges displayed in Table 1. Note that thickness step was set to 2.825 Å in all simulations, and both values in each range in Table 1 were included. We labelled the simulations from 0 to 85535, thus it is possible to obtain the layer sequence directly from our nomenclature and Table 1, more info is provided in the Appendix. Figure 3 shows the material gain for all simulations performed with x = 0.24. All structures that had material gain >15cm −1 were selected into a subset over which we conducted additional simulations by restricting number of states per period to investigate the leakage effects. Figure 3 does not provide much information on design quality, however, since the layer variation was constructed through nested loops, we can observe the effect of variation of each layer in Table 1. For example, the first layer was varied in nine steps and every section depicted by red dotted lines in Fig. 3 corresponds to a different thickness in the range 50.85 -73.45 Å. If we enlarge one of these sections, we would get 6 regions that correspond to variation of the second layer. The choice of ranges in Table 1 was made so that we observe a parabolic envelope every time we enlarge the layer effect in order to ensure the correct search for the optimal design.  Table 1. Table 2 presents the designs found to be promising candidates. We conducted a very large amount of simulations, and the best 10 designs for each Al molar fraction are available in Supplement 1, along with more detailed information (best designs per each molar fraction).
We will use the data with most states per QCL period in discussion, however we do acknowledge that the highest state may have ended up in the continuum, making some data unreliable.
One of the observations made in simulations was the appearance of designs with very high material gain, while there was no transition to support the frequency corresponding to that gain (offsets over 1 THz). These simulations were filtered as false results and we set the tolerance of 200 GHz. This tolerance was set only in the initial simulation with four states and, as Table 2 shows, some designs have offset over 200 GHz. This offset means that DM model reports the lasing frequency to be f , however the lasing transition should normally correspond close to the energy difference E 32 .
In our experience with the DM model, this is a somewhat common effect. Similar behaviour occurred in the former record temperature design of 200 K [12] and was attributed to competing oscillatory strengths and the fact that LLL was formed of two narrowly spaced states. In our simulations of 200 K structure, the DM model yields a frequency of 3.13 THz and energy separation of 3.3 THz at the peak of material gain dependence on bias, while the experimentally observed frequency was 3.22 THz (at threshold), thus 3.3 THz is a more reliable result due to Stark effect. However our simulation of 210 K structure [13] yielded a frequency of 3.88 THz and energy separation corresponded to 4.15 THz, while the structure lased at around 3.88 THz at peak power [13] (which most likely corresponds to the peak of the material gain in our simulation).  This offset is caused by the other states present in the design and, interestingly, is proportional to the anticrossing energy between resonant states that DM model does take into account, and we will rely on using the frequency obtained by DM model, with caution towards very high offsets from the transition frequency corresponding to E 32 .
Our first criteria was finding a design with the highest value of material gain around 3.5 THz, where we expect the minimal loss. Our second criteria was having a design with low electrical heating which we roughly described by H F introduced in Eq. (2), simply as a product of current density and resonant bias. To join these into a single target, we initially introduced a quality factor Q simply as the ratio of material gain and H F . However, in the simulations we have noticed somewhat unusual behavior in current density dependence. In our experience with QCL modelling, a change of monotonicity in current density dependence on electrical bias, commonly referred as non-differential-resistance (NDR) region, typically causes a sudden break in lasing performance [13], even though theoretical models (like our) may still indicate the positive net gain which was the case in some of our simulations. The sudden termination of lasing in NDR has been attributed to the fact that, in typical LO-phonon structure, the tunnelling process stops at some bias and since QCLs are mainly driven by the current, a sudden shift in voltage occurs in order to find a better "channel" to conduct the driving current, therefore breaking the lasing process. The design in Fig. 2 has two resonant tunnelling processes, and it is possible that it may be more resilient to sudden breakdown of lasing. If for instance ILL 2 → ULL stops, but LLL → ILL 1 is still ongoing, the entire design behaves as the scattering assisted structure presented in Fig. 1(c)). This property of our design needs experimental verification, which is why in Table 2 we present the model output at two bias points: at the peak of material gain dependence (K) and at the peak of the current density dependence (which occurs at bias K n ). This issue is negligible in Table 2, but will be further discussed in four-well design analysis.
The data in Table 4 leads to several conclusions: • The inclusion of higher states in the simulation decreases the material gain and increases the current density. This illustrates the leakage processes due to higher undesired states present in the simulation.
• The previous effect is significantly weaker in designs with higher barriers, as expected.
• The Q factor is somewhat higher for higher barriers. This was not initially expected, as higher barrier designs should have lower material gain, since the sheet doping was constant in all simulations. The main reason that somewhat opposite effect occurred was most likely due to thinner injection barriers that were used for x = 0.23 − 0.25.
• The lasing frequency is more stable with inclusion of higher states in simulation in designs with higher barriers, and there is a significantly larger number of designs around 3.5 THz (we found none with high material gain in simulations with x = 0.15).
• The structures whose peak material gain does not occur at the peak of the current density vs. bias dependence are present, but not significant in designs of our interest (around 3.5 THz).
The discussion of results in Table 4 is nearly identical to the discussion of two-well LO-phonon THz QCLs [13]. The main difference however, is in the fact that we used sheet doping density of only 2 · 10 10 cm −2 , smaller than in the record structure [13] that uses 4.5 · 10 10 cm −2 , and we can improve the material gain by either increasing the doping or using a thinner injection barrier. We found that x = 0.23 − 0.25 provides the best designs since the leakage due to high states is negligible, the NDR gain difference is also negligible, and Q factor is the highest. We will now further analyze these designs: The structures presented in Fig. 4 have very similar layer sequences and they all share similar properties in terms of transport mechanisms. The first four states are colored as in Fig. 2, to reflect that the goal for finding such design has been achieved, although ILL 1 is localized more centrally for x = 0.24 and x = 0.25 in the period. As discussed in Section 2, the simulation yielded best designs with LO-phonon separation above 36 meV, which was also one of the main design arguments in [13] that prevents LO-phonon leakage. Interestingly, all designs display a very interesting formation of higher, undesired states in the period. Both in D_3_24_81341 and D_3_25_79626, the fifth and the seventh state are spaced by approximately LO-phonon energy, meaning that these two undesirable states would additionally assist the ILL 1 state, thus offering a different paradigm to the design in Fig. 2. This is present to an extent in the lower barrier design D_3_16_81341 as well. Note that the seventh state in this design is critically close to the continuum, and Table 1 shows a distinct deterioration of the material gain when this state is included in the simulation. In our experience, including a continuum state in the simulation often causes the Schrödinger-Poisson solver to break down, because some matrix elements needed for LO-phonon scattering calculation would have non-physical values, however, a state critically close to the continuum allows the simulator to calculate the transport. We cannot asses whether the sudden drop in material gain reflects an actual physical effect, or is a numerical issue. It is possible that only six state simulations should be considered, however, in our experience, including higher states in simulation of resonant-phonon structures typically does not deteriorate the material gain as significantly as in Table 1 and we decided to keep the analysis with seven states per period, as it is a worse scenario.
The best structure in three-well simulation is D_3_25_79626 (as it was also promising in x = 0.24 simulations), however we would recommend for D_3_16_81344 to be experimentally tested as well, since the six state simulation data are promising, and x = 0.16 structures may have smaller growth fluctuations than x = 0.23 − 0.25 structures. We also recommend changes in doping (and injection barrier widths), whose effect is illustrated in Fig. 5. The model shows a linear increase of current density with doping, the material gain on the other hand undergoes the saturation effect, while frequency is not significantly affected. Since a high current density would cause higher electrical heating, the optimal value needs to be found experimentally. The saturation effect of material gain in Fig. 5 suggests that value in the range 2.5 − 3 · 10 10 cm −2 is potentially optimal.

Four-well design
The simulations with x = 0.18 − 0.23 in three-well structures did not provide significantly promising designs. This can be understood through a closer look at states five and six in Fig. 4. For lower barriers, ILL 1 is constructed properly, however the fifth state is in close resonance with ULL from the previous period, thus creating an inefficient dual lasing channel. The fifth state rises with the increase of barrier height, where the optimal distance from ILL 1 would be the LO-phonon energy, which is the reason why designs with x = 0.24 and x = 0.25 have shown the best performance.
The previous temperature record of 200 K [12] used a three-well structure with x = 0.15. This was a resonant phonon design, where there was another level below LLL that was assisting the  extraction process. It is clear that detrimental effect of the fifth state for lower barriers can act beneficially if we implement the paradigm in Fig. 2 in a four-well structure (simply find a layer sequence that adds ILL 1 state in typical three-well LO-phonon design). Another major advantage of such four-well structure would also be in a significant reduction of resonant bias and current density.
We therefore conducted additional set of simulations of four-well designs by altering all layers' widths (apart from the injection barrier) and keeping the sheet doping density to 2 · 10 10 cm −2 , as shown in Table 3.  36.725, 33.9 Å, respectively. Only the widest well has been doped. 73.45 Å of the well was undoped, while the rest of the well's width was doped so to correspond to 2 · 10 10 cm −2 sheet doping density for all Al mole fractions. The lattice temperature was 250 K in all simulations.
In Table 4 we present the designs that were found as promising candidates with similar criteria as discussed for three-well designs. The 3.5 THz designs were our main focus, which are not necessarily the designs with best Q factor or material gain. In Supplement 1, Dataset 1 [41], Dataset 2 [42], Dataset 3 [43], and Dataset 4 [44], we provide a more extensive simulation tables and simulation results. The only difference in the simulation procedure is that now we are looking for structures with five states, where LLL state would have an additional assisting level for efficient extraction. This means that ILL 1 would be the fifth state in the period, we can still view this design with paradigm in Fig. 2 where we allow LLL to be effectively a miniband instead of single quasibound state. To account for leakage effect, we first made simulations with five states per period, selected promising candidates and then conducted additional simulations on that subset by allowing six, seven and eight states per period. Table 4 typically displays higher Q values than Table 4. This is mainly because of lower values of resonant bias and current density than in three-well structures. The table values offer similar conclusions as discussed previously, the difference is that in some designs the material gain has sudden drop when the sixth state is included in the simulation, and some designs may have slightly higher gain with eight states in a period. This occurs due to more complicated structure where the undesired states (six, seven, eight) interact with ULL or with both states that form the LLL. We also found that most designs in Table 4 have the sixth state acting as ILL 1 , instead of the fifth state.
The issue with the frequency offset between DM prediction and E 43 energy is also present, however cases with very high offset may indicate the lasing states do not correspond to E 43 difference. Similarly, as in three-well simulations, frequencies around 2.6 THz are more favoured    Table 4 share similar behavior and excellent performance, making the selection of the optimal structure challenging. Our main aim is low current density, high material gain and low resonant bias. The success of these designs lies in the fact that parasitic state which is in resonance with ULL from the previous period has virtually no effect on transport. The three-well designs needed higher barriers to "shift" this state away from ULL (Fig. 4). An interesting deviation from paradigm in Fig. 2 is that ILL 1 is the sixth state in period and it is separated from the parasitic fifth state by a similar energy as that between the lasing states, thus creating a potential dual lasing channel, and the energy difference between ILL 1 and ULL is significantly higher than energy of LO-phonon resonance and of the order of ∼ 50 meV. The paradigm in Fig. 2 is still satisfied, as LO-phonon process may occur on that energy scale at very high temperatures [30], however it is also possible that the structure operates with two lasing processes, although the upper transition is highly diagonal and it is somewhat unclear how the fifth state is depopulated even though it is ∼ 36 meV from ULL.  Figure 7 shows that all designs we presented may slightly benefit from thinner injection barriers than the ones we used. The highest material gain could be achieved with 11 monolayers in D_4_0.23_103899 instead. However, designs with thinner injection barriers have high material gain in NDR region and only with 13 and 14 monolayer thickness (and more) the peaks of the current density and material gain are nearly aligned. Discussing this effect further is outside of the scope of this paper. We used typical injection barrier thicknesses and did not vary the doping level in our simulations, mainly because both parameters are obtained through somewhat empirical approaches and our model is tight-binding which may have placed some unrealistic constraints in determination of the electronic structure. Note that we also did not use any preset value of loss, thus we did not calculate the optical power and the corresponding increase in current density dependence when the structure is lasing. In our experience with the DM model, QCLs normally do not provide material gain in NDR region and the current density dependence is somewhat imprecise [37]. In the simulation for x = 0.23 we used injection barrier width of 12 monolayers, and using thicker barriers may also provide lower values of current density at the cost of some loss of material gain. For instance, with 15 monolayers thickness the current density drops from 1270 to 1019 A cm 2 at the cost of material gain dropping from 14.7 to 11.9 cm −1 , which is not too detrimental because this could be rectified with higher doping in order to achieve a structure where high gain is not located in NDR region. However, this effect needs to be verified experimentally, as it is possible that our design may operate different from other THz QCLs since it undergoes two resonant tunnelling processes.

Cut-off temperature performance
The five designs we presented so far all satisfy the design criteria in Fig. 2, however the simulations were performed at 250 K, generating the structures that at first sight do not appear convincingly better than the two most recent record temperature structures.
The 200 K structure has the loss evaluated to be around 20 cm −1 , while the 210 K structure was designed for a higher frequency and authors [38] predicted the loss of 30 cm −1 , making the DM prediction very precise in Fig. 8(a). However the latter prediction may have been too strict, as this structure displays material gain above 20 cm −1 nearly up to 240 K, which is in line with our design predictions as well in Fig. 8(a). We have two main arguments that should be considered in attempting our designs: • The sheet doping density in the 200 K structure [12] was 3 · 10 10 cm −2 , while in the 210 K structure [13] it was 4.5 · 10 10 cm −2 . In all the designs presented in this work, the sheet doping density was set to 2 · 10 10 cm −2 , which is significantly lower and indicates that all our designs would provide a higher material gain than the current record structure. Figure 8(b) clearly illustrates this. All the designs in Fig. 8(a) have the same sheet doping density as 210 K structure. Note that this has significantly increased the current density of our designs in Table 2 and 4, roughly by a factor of 2.25, and N s = 4.5 · 10 10 cm −2 may not be the optimal doping value due to gain saturation effect displayed in Fig. 5. We only did this so we can clearly illustrate the quality of our designs over the 210 K structure.
• Figure 8 also displays the slope of gain deterioration, where we can notice that three-well designs have slopes between those for the two referent realized structures, and mainly provide high gain due to a high gain at 10 K. The four-well structures have significantly better slope than the current record structures, and a great argument can be made for growing D_4_0.18_79421. (a) Our structures have the sheet doping density of = 2 · 10 10 cm −2 , while 200 K structure has = 3 · 10 10 cm −2 and 210 K structure has = 4.5 · 10 10 cm −2 .
(b) All structures have the sheet doping density of = 4.5 · 10 10 cm −2 . Fig. 8. The temperature dependence of peak material gain for two recent record structures of 200 K [12] and 210 K [13] and several designs from this work. Each simulation was conducted at respective resonant bias for each design. Interestingly, the four-well designs with higher barriers provide the best slope of gain deterioration, however designs presented in Fig. 8 do not show very high gain at low temperatures. This can be attributed to diagonality of lasing transition and energy separation between states that undergo LO -phonon transition. Design D_4_23_103899 has E 21 =48.7 meV, while D_4_18_79421 has E 21 =43.3 meV. Since the LO -phonon scattering mechanism is not dominant at low temperatures, design D_4_18_79421 would have a better optical power performance than D_4_23_103899 at low temperatures, as displayed in Fig. 8, however as LO -phonon interaction is dominant at high temperatures, design D_4_23_103899 is becoming more robust. For that reason we presented the performance of additional designs in Fig. 8 as trade-off solutions, making D_4_22_129641 potentially the best candidate, as it behaves nearly identically as D_4_18_79421, with slight benefits in parameter values in Table 4. Note that designs with high LO-phonon energy separation may still have high gain at low temperature if the lasing transition is more vertical, which is the case with D_4_22_129641 which has E 21 =45.6 meV.

Conclusion
In this work we have reviewed the existing THz QCL designs by creating an effective three-and four-level schematics that underpin the key mechanisms and differences between experimentally realized structures. This review led us to a novel idea for a four-level structure that undergoes two resonant tunnelling process and hybridizes scattering-assisted and resonant-phonon QCL designs into a design that may take advantage of commonly argued processes for temperature degradation, so to act beneficially in improving the lasing performance.
The degradation of material gain at high temperatures in resonant phonon structures was argued mainly to occur due to the fact that LO-phonon process may depopulate the ULL as well, since this process is active far after the resonance energy difference at 36 meV (in GaAs). Another process that was often argued as the gain degradation mechanism was the interaction of lasing states with undesired higher states from adjacent period (i.e. "thermal backfilling"). Both these issues had a solution in using higher barriers, at the cost of lower material gain, that can be compensated by a higher doping density at the cost of higher current density, which however leads to detrimental electrical heating of the device. The current record temperature design [13] addressed all these issues by finding the compromising trade-off solution.
The design proposed in this paper in Fig. 2 attempts to set the higher states in QCL period, that are involved in thermal backfilling, to be in LO-phonon resonance with ULL below. This creates a symmetric structure in terms of transport processes that provide more efficient population inversion maintenance by collecting any potential leakage effects into another "tunnelling + LO-phonon transition" cycle.
Majority of THz QCLs were designed by similar methods as MIR structures, where one could design the active region (lasing states) and collector/injector (typically ILL) separately. This comes at some cost of neglecting the 'compound' effects when these sections are joined. Additionally, the design procedures typically set the resonant bias beforehand and seek structures with particular energy differences, and evaluate the design quality through oscillator strengths which are proportional to the product of energy difference and square of dipole matrix element (which is proportional to the wavefunction overlap). This procedure has led to lengthy discussions on whether higher (vertical transition) or lower (diagonal transition) dipole elements between the lasing states provide a better trade-off solution. Such techniques were historically necessary due to models with only a few effective states, that were inherited from MIR QCL design approaches, and only recently [38] a more powerful transport models such as NEGF are being used in designs which yielded the current record temperature performance [13].
The NEGF model does provide high quality information on electronic structure, however it suffers from very high numerical cost. In the numerical section of this work, in our opinion, we performed one of the most extensive analyses of THz QCLs, by using significantly simpler DM model. Our model is the generalization of the DM approach used to optimize the former 200 K [12] record design, and we observed a high quality predictions of cut-off temperature that was also reported by authors that used a very similar model [37]. Due to very low numerical cost of our model, and the access to supercomputer cluster, we were able to simulate a very large number of THz QCL structures and perform brute-force technique in seeking for the optimal structure.
We do acknowledge that an optimization algorithm would have been a better approach, however, similarly as in [12], the DM model cannot account for empirical effects and design imperfections, electric field domain formation, or fluctuations of layer thicknesses. Additionally, the results in Fig. 8 do not visually show the full extent of quality of proposed designs over the 210 K record structure as the doping density is kept at very low level in our simulations. We led the discussion in this paper that electrical heating needs to be rectified, which may not necessarily be the requirement for high temperature operation, as the 210 K structure [13] is very highly doped and has very high threshold. Therefor each design that we presented may provide significantly more enhanced characteristics if electrical heating is ignored by adopting that duty cycle may be engineered to be very low.
The main aim of this work is to present a novel design idea and we hope that similar findings could be verified by more extensive models, such as NEGF, that would support attempting this design experimentally. Our main argument is that LO phonon transitions should not be design around 36 meV if high temperature performance is needed and that a very careful consideration of states above ULL needs to be considered. The best designs in our simulations are not just the one that placed state above ULL into LO phonon resonance, but also placed the next level above at LO phonon resonance ( Fig. 6(b)) meaning that ideal design should form all levels above ULL in LO phonon resonance in order to pump ULL most efficiently. Note that such designs also show a smaller parasitic leakage because the subsequent LO phonon transitions above ULL are contributing to pumping of ULL and not to leakage processes.
We also used potentially disputable arguments for focusing on the designs at 3.5 THz, that may not have the lowest loss or defining the Q factor as the ratio of material gain and the product of current density and resonant bias instead of setting a predetermined value of loss and using net material gain value instead. For this reason, we have also presented other prospective designs that our simulations yielded, and provided an extensive data set from our simulations in Supplement 1. This data contains high level of detail and offers extended set of parameters to those presented in Table 1 and 3 and may be used in seeking devices with different frequency, current density, resonant bias and material gain (at 250 K).
first layer corresponds to the 6th point, we would need to find a whole division of (65320-6*9504)/1584=5.23, which corresponds to the 5th (the last) point in the second layer range which is 19.775 Å. The next layer was simulated in 12 points, in the range 19.775 -50.85 Å, thus the next layer in the sequence is (65320-6*9504-5*1584)/(1584/12) = 2.84 which corresponds to 25.425 Å, the next layer is in the range 5.65 -33.9 Å of 11 points and by similar analogy (65320-6*9504-5*1584-2*132)/12 = 9.33 which is 31.075 Å and the final layer is simply (65320-6*9504-5*1584-2*132)%12 = 4 which is 144.075 Å. Therefore, the design D_3_0.15_65320 has the layer sequence 67.8/19.775/25.425/31.075/144.075 Å (where layers in bold are barriers), the injection barrier for the 0.15 design was 42.375 Å thus this sequence needs 21.1875 Å at its ends to fully describe the layers of this design (and, as we mentioned, the widest well was doped, 84.75 Å of it is not doped, while the rest is doped to correspond to 2 · 10 10 cm −1 sheet doping density). The second approach in decoding the layer sequence from this nomenclature would be to notice that the terms in this calculation as (65320-6*9504-5*1584) are simply modulo division of (65320-6*9504)%1584 and the decoding procedure can be easily numerically implemented.

A.2. Thermal constant
The ridge of QCL structure is typically grown on top of a ∼200 µm thick substrate where the cold finger temperature is set on the bottom contact (at the substrate side). The electric power is dissipating mainly in the ridge, thus the temperature is expectedly rising towards the top contact (on the top of the ridge). Since the transport is occurring in one direction, the first approximation for the heating of the device is given by the Fourier law of conduction [30]: where ρ m is the material density, c p is the specific heat capacity, k z is thermal conductivity in lateral direction, P is the electrical power and V QCL is the volume of the active region (we assume that power is only dissipated in the volume of the ridge). In steady state, this equation only requires the information on thermal conductivity which is spatially dependent due to the heterostucture layers in the active region. However, for illustration let us assume that both k z and P are constant and that bottom of the ridge is at same temperature as the cold finger (i.e there was no heating through the substrate) and set this point as z = 0, so ∂T ∂z | z=0 = T 0 and that the top of the ridge (z = H QCL ) is in the vacuum T | z=H QCL = T sub . The analytical solution )︂ . We can now obtain effective active region temperature by integrating over z and dividing by ridge height, which results in T = T 0 + H QCL P 3k z WL c (note that we used V QCL = H QCL L c W where L c is the cavity length and W the substrate thickness).
An empiric approximation for heating effects in QCL structures of the form: T = T 0 + R TH P where R TH is a thermal constant, is the usual model for thermal effects in QCLs [30] and it is clear that the R TH constant is related to volume of the QCL ridge (as it originates from term P V QCL in Eq. (3), in general).
The analytically derived thermal resistance under approximation (and averaging) we made, shows that R TH = H QCL 3k z WL c and this gives somewhat intuitive conclusion that QCLs with wide ridge, long cavity and shorter active region display better thermal performance. However, the electrical power is the product of current and voltage, which when re-scaled to the current density and electrical field P = IV = JKWL c H QCL cancel the dimension terms in R TH and only suggest that the height of QCL ridge H QCL has the main effect on heating (however this conclusion holds if distribution of electric bias and current density is uniform). Although our analytical derivation has been obtained through ambiguous approximations, it may be generally concluded that in Eq. (2), the product of R TH and ridge volume is somewhat canceling the effect of the ridge volume as discussed in Section 2.