Doping study of two-well resonant-phonon terahertz quantum cascade lasers part I: doping profile dependence

The influence of impurity doping on GaAs-based two-well resonant-phonon terahertz quantum cascade lasers is investigated theoretically, and efficient doping schemes are discussed. By using the rate equation model, the impacts of dopant amount, position, and distribution on the performance of a high-performance device is simulated focusing on a single module. The calculated optical gain is found to have a peak over the range of sheet doping density from 1.0 × 1010 to 1.0 × 1012 cm−2 in all eight doping conditions examined in this work. Among these patterns, the devices with the undoped condition and homogeneous-doping in phonon-wells mark high optical gain, and the latter is also resistant against the detuning of subband alignment due to band-bending under the high doping conditions. Furthermore, based on the simulation results, a modulation doping scheme whose active cores include both doped and undoped modules is suggested and discussed.


Introduction
Quantum cascade lasers (QCLs), first proposed by Kazarinov and Suris in 1971, 1) have enabled lasing in mid-infrared (MIR) 2) and terahertz (THz) 3) frequency regions and have attracted a great deal of attention because of their potential to realize various state-of-the-art sensing systems and advance wireless communication technologies. [4][5][6] Compared to conventional diode lasers, 7,8) which emit light by the recombination of electrons from a conduction band and holes from valence bands, QCLs have three distinctive features: wavelength tunability, sharp lasing spectra, and carrier recycling. 1,9) These characteristics are attributed to fundamental QCLs' driving mechanisms by which carrier transport for radiation occurs by using only electrons. First, diode lasers and QCLs commonly create population inversion between two bands or subbands for stimulated emission, and light whose energy corresponds to the energy difference between the two bands is generated. Therefore, the wavelength of light for diode lasers mostly depends on the bandgap, which is unique to materials and has little flexibility for tuning. On the other hand, the wavelength of QCLs is more flexible in design through the adjustment of layer thickness and/or barrier height (Al-composition in the GaAs/AlGaAs material system). Secondly, the optical linewidth is dependent on the curvatures of the energy-momentum parabolic curves (the E-k dispersion curve). The optical transition of diode lasers happens between a concave conduction band and convex valence bands where carriers take a certain range of momentum, resulting in optical linewidth broadening. In QCLs, on the other hand, the electrons contributing to optical transitions reside only in the subbands inside a conduction band whose band-curvatures are almost the same. Thus, the optical linewidth of QCLs is narrower. Finally, in QCLs, it is also noteworthy that electrons are not annihilated by radiation and are instead recycled for the next optical transition. This carrier recycling effect is beneficial in enhancing internal quantum efficiency.
Thanks to these superior features, MIR QCLs, first demonstrated by Faist and Capasso et al. in 1994, 2) have been commercialized for a wide range of applications since the realization of their room-temperature operation. 10,11) On the other hand, THz QCLs were demonstrated by Köhler et al. in 2002, 3) but their highest operation temperature remains 250 K even at present and a solution to device performance degradation at high temperatures is urgently required. 12) To address this issue, the design of active cores and waveguides is important. As for waveguide designs, metalmetal waveguides, especially Cu-Cu waveguides, are confirmed to be the most effective for low threshold gain and high operating temperatures. 13,14) However, we are still seeking appropriate design guidelines for active cores that will account for a variety of features: carrier transport schemes, 3,12,[15][16][17][18][19][20] barrier height, [21][22][23][24][25] structural complexity, [26][27][28][29][30] material systems, [31][32][33][34][35] and impurity doping.  Among these components, a clear design guideline for carrier transport scheme is found by Bosco et al.: active cores composed of fewer lasers are advantageous for high-temperature operation. 15) Actually, the current highest operation temperature is demonstrated by a device with only two wells. 3) On the other hand, design guideline for other components is still unclear. Thus, this study focuses on the influence of doping methods on device performance.
Doping has been studied for both MIR-and THz-QCLs  with seven factors: ionized-impurity scattering, [43][44][45][46][47][48] doping amount, [36][37][38][49][50][51][52][53][54][55] doping position, 18,39,40,57,58) doping distribution, 59,60) device polarity, 41,61,62) band-bending, 42,63) and free-carrier absorption. [64][65][66] The impact of ionized-impurity scattering on carrier transition rate and optical linewidth is the most important feature in doping studies. [43][44][45][46][47][48] In terms of carrier transition rate, ionized-impurity scattering is not a dominant phenomenon in THz QCLs compared to LO-phonon and interfaceroughness scattering. However, in accordance with Fermi's golden rule, where the scattering rate is determined by an overlap amongst wavefunctions of the initial and final states and dopants, when the doping position is close to the lasing sites, non-radiative scattering rate can increase due to ionized-impurity scattering. 39) This effect is less influential in MIR QCLs and early-stage THz QCLs. Because a doping region (injection region) and an active region are spatially separated in these devices, the carrier transition between lasing subbands is less likely to be exacerbated by ionizedimpurity scattering. In recent THz QCLs composed of fewer wells, on the other hand, donor impurities are doped very close to lasing sites. In this case, the degradation of population inversion is inevitable. This issue is thus thought to be more critical in short two-well THz QCLs. 12,15,19,30) As for optical linewidths, ionized-impurity scattering is one of the main causes of optical linewidth broadening. 47,48) Therefore, ionized-doping greatly influences optical gain through optical linewidths.
Studies on an optimal doping amount have been the most frequently reported, and representative works for THz QCLs are listed in Table I. In these previous works, 47,[49][50][51][52][53][54][55][56]63) resonant-phonon (RP), bound-to-continuum (B-to-C), and split-well direct-phonon (SW-DP) structures are used with a different range of sheet doping densities from 4.3 × 10 9 to 3.5 × 10 11 cm −2 , and the highest operation temperature, optical gain, and output power are reported to have optimized values. 49,50,55) This fact is interpreted to mean that dopants are necessary to obtain a sufficient optical gain to surpass the threshold gain; however, too high a dopant injection causes high free-carrier absorption [64][65][66] and band-bending, 42,63) resulting in increased waveguide loss and undesirable misalignment of subbands, respectively. Free-carrier absorption must be suppressed in THz QCLs because it is proportional to the square of dopant density (N d 2 ) and wavelength to the power of three-and-half (l 3.5 ). 64,65) Despite relatively few studies, the influence of doping position has also been investigated. 18,39,57,58) Previous research using MIR QCLs reveals that when dopants are injected in the vicinity of the upper lasing state (ULS), device performance degrades remarkably. 39) In THz QCL research, Wang et al. also examined the influence of doping positions on the highest operation temperatures, using threewell resonant-phonon devices. 58) In comparison between a phonon-well doped device and a radiation-barrier doped one with a typical sheet doping density (3.0 × 10 10 cm −2 ), the demonstrated highest operation temperatures were almost the same. Demić et al. also theoretically confirmed that device performance is independent of doping position with 2.0 × 10 10 cm −2 of sheet doping density. 18) On the other hand, Grange, through his theoretical work, demonstrates the possibility of improving device performance with different doping positions. 57) Because this type of research is expected to be influenced by other components-doping amounts or distribution-more-detailed research is necessary.
Lastly, a few studies on doping distribution have also been reported. The most frequently used doping method for recent resonant-phonon devices is narrow-rectangular doping in the center of phonon-wells. 12,15,16) On the other hand, δ-doping method, which can minimize interactions between doping regions and wavefunctions, is also thought to be effective at enhancing device performance. 59,60) In experiments, furthermore, the influence of dopant migration should also be considered. During epitaxial growth, dopants diffuse into more freshly grown layers. If these diffused dopants reach radiation areas, device performance is degraded. 61,62) To avoid this issue, the order in which layers are grown should be carefully considered.
This study aims to theoretically examine doping methods for two-well resonant-phonon THz QCLs and to provide doping design guidelines for high-temperature operation. First, based on the doping-related phenomena described in Sect. 2, the optical gain of eight doping patterns is calculated for a single module by using device simulation, and the best doping pattern is determined and analyzed in Sect. 3. In Sect. 4, an effective doping scheme is suggested and discussed, including phenomena in waveguides that are not dealt within the simulation.

Model
The device simulator based on the rate equation has been developed in our previous work. 25,67) In calculating carrier scattering rates, LO-phonon (LO), interface-roughness (IFR), ionized-impurity (IMP), alloy-disorder (AD), and electronelectron (EE) scatterings are included with a screening effect model. For the optical linewidths inside a module and pure-

032002-2
dephasing times between modules, the four scattering phenomena-except for electron-electron scattering-are included. The carrier transport between modules is calculated by tight-binding theory. To calculate electrical current and effective subband carrier density while taking care of carrier leakage from bound states to continuum states, third-order tunneling current theory is also applied. 25) As device parameters directly influenced by impurity doping, current density, optical gain, carrier scattering rate, optical linewidth, and band-bending have often been discussed. In addition to these direct influences, indirect effect on some of these parameters caused by band-bending is also important to consider simulation and experimental results. On the one hand, band-bending issues subband detuning at lasing alignment conditions. 63) Such potential distortion, on the other hand, affects other device parameters such as oscillator strength and tight-binding parameters via wavefunction and eigen energy calculated by Schrödinger equation, resulting in errors in current density and optical gain (Fig. 1). To understand such complicated influence of impurity doping seen in the later sections, the theoretical models of scattering rate, optical linewidth, dopant activation ratio, and the δdoping distribution model are introduced in this section. To simulate band-bending, Hartree approximation is used.

Ionized-impurity scattering models
Ionized-impurity scattering is the most important physical phenomenon in this doping study, and the carrier scattering rate and optical linewidth due to it are calculated by Fermi's golden rule 68) and Ando's theory, 44,69) respectively, using the same scattering Hamiltonian˜( In Eq. (1), q 0 and e ( ) z s denote the elementary charge and the static permittivity, respectively. The locations of a scattered electron and ionized dopant are also described by xy,0 0 When electrons transfer from an initial subband ñ |i to a final subband ñ | f , the scattering rate is calculated by Fermi's golden rule described in Eq. (2) In Eq. (2), electrons' momentum and total energy are described by k i and k f and by E i T and E , f T respectively. The reduced Prank constant is described as . Based on Eqs. (1) and (2), the ionized-impurity scattering rate is calculated as in Eq. (3), and thermally averaged rates are used in the rate equation The effective mass of electrons is denoted by m , * and Q W and J W are respectively calculated by Eqs. (4) and (5) Next, Ando's theory, derived from the first principal calculation, is described in Eq. (6). In it, the optical linewidth of a transition between ñ |i and ñ | f is denoted as G . if Instead of the total energy, kinetic energies E i and E f are used in the delta function. By a similar mathematical treatment, the optical linewidth due to ionized-impurity scattering is obtained as in Eq. (7). Like the calculation for scattering rate, components G Q and G J are calculated by Eqs. (8) and (9). The difference between Eqs. (4) and (8) comes from the differences in conserved energy types For calculating the Coulomb scatterings, the influence of screening effect is also indispensable, and it is included by replacing the relative wave vectors Q W and G Q with the ones that include the screening wave vector K sc as in Eq. (10) 67) Even though the optical linewidth is calculated by using Eq. (7), the behavior of this parameter is a little difficult for us to imagine from the theoretical form. To obtain a crude sense of the relation between optical linewidth and other parameters such as wavefunction and doping, a series of example calculations is implemented based on a GaAs/Al x Ga 1−x As-based double-quantum-well (70/20/140 in Å. Bold is the barrier. x denotes the Al-composition), and the results are summarized in Fig. 2 where ñ |i and ñ | f denote the ground states in each well. Calculations are implemented with a sheet doping density 3.0 × 10 10 cm −2 and population fractions (r , i r f ) = (0.9, 0.1) under no electrical field and 200 K of heat-sink temperature conditions. Unlike the full rate equation model, some physical phenomena such as the Hartree potential and electron temperature model are ignored. In calculation, the wavefunctions of the initial and final states are fixed, and the 20Å-rectangular doping region is shifted among four different positions described in Fig. 2(a).
In the results, the scattering rate and optical linewidth take the highest values in Position C. In this condition, the wavefunction of the initial state forms a peak in the left well, and it has a marked overlap with the doping area. The scattering rate decreases as the doping position gets further from the left well. In Position B, the optical linewidth is one order smaller than the ones of other conditions. This is because the optical linewidth is mainly determined by an overlap of the doping area and the probability density difference of the two wavefunctions as described in Eq. (9). Furthermore, the same tendency is observed even when population fractions and barrier height (Al-composition) are changed as in (r , i r f ) = (0.5, 0.5) and (0.2, 0.8), and the Alcomposition changes from 15% to 30%.

Dopant activation model
As is often the case with QCL simulation, carrier density is assumed to be the same as doping density. However, even at room-temperature, not all dopants are activated. So, a selfconsistent calculation is necessary to estimate the net carrier density. In the case of n-type semiconductor materials, the activated dopant density + N d depends on the donor level E D generated between the conduction band E C and Fermi level E F due to impurity doping and is described by Eq. (11) In Eq. (11), the three-dimensional doping density is denoted by N , d and k B is the Boltzmann constant. As seen in this equation, the dopant activation ratio is influenced by the lattice temperature T .
L For the band diagram in Fig. 3, the energy difference between the Fermi level and donor level is described by Eq.
The ionization energy of doped Si in GaAs is described by -E E , C D and its value is 5.85 meV. 70,71) -E E C F is the energy difference between the conduction band and Fermi level and is obtained by Eq. (13) by using the effective density of state N C in Eq. (14) Finally, the value of the net three-dimensional carrier density n 3D is assumed to be the same as the calculated activated doping density described in Eq.  Having said that, to calculate + N d by Eq. (15), n D 3 is required in Eq. (13). Therefore, these calculations are selfconsistently processed with an initial assumption of » n N .

D d 3
The dopant activation ratios of Si in GaAs at 200 K and 300 K are approximately 84% and 91%, respectively. 71)

Dopant distribution model
In this study, conventional rectangular doping and δ-dope distributions are examined considering the back-ground doping (1.0 × 10 14 cm −3 ), and these doping distributions are determined by simple theoretical models and the sheet doping density N . In rectangular doping, dopants are doped homogeneously over a range of L dope described in Fig. 4(a), and the sheet doping density is calculated using a local doping density N .
In actual situations, segregation happens as depicted by a wiggly line in Fig. 4 For δ-doping, an experimental model considering dopant migration is adopted as described in Eq. (18). 59,60) In Figure 4 is a peak of doping density and z p is the position of the doping peak. Dumping coefficients of exponential curves are denoted by a and b. In this study, a = 4 and b = 24 are extracted from the reference data. 60) From Eq. (18), the sheet doping density is calculated as in Eq. (19).

Simulation
Based on one of the latest two-well resonant-phonon structures, the influence of doping conditions on device performance is computationally investigated. The base structure (G652) has been designed by Khalatpour et al. 12) using a GaAs/Al 0.30 Ga 0.70 As material system (the top figure in Fig. 5), and the carrier transition for lasing happens amongst the lowest three subbands ( ñ |1 ∼ ñ |3 ). At a lasing alignment condition, the injection state ¢ñ |1 aligns with the upper lasing state (ULS) ñ |3 to inject electrons to the next module, and population inversion is created between ULS and the lowerlasing state (LLS) ñ |2 . The electrons in LLS are extracted by fast LO-phonon scattering from LLS to the injection state. The basic performance of this structure has already been investigated in our previous work. 25) Doping conditions examined in this study are divided into eight patterns depending on positions and distributions (the bottom figures in Fig. 5). To investigate the influence of doping positions, narrow-rectangular doping patterns in phonon-wells (A), radiation-barriers (B), injection-wells (C), and injection-barriers (D) are simulated. Following that, δ-doping in phonon-wells (E), homogeneous-doping in phonon-wells (F), and all layer doping (G) are examined to observe the influence of doping distributions. Furthermore, the undoped condition (H), where electrons are assumed to be injected from an external module, is also simulated. In this condition, the equations in Sect. 2.3 are not used. Sheet doping density is constant (2.69 × 10 8 cm −2 ), and (injected) sheet carrier density is swept as a variable. In any conditions, n-type natural doping is assumed to be 1.0 × 10 14 cm −3 over all regions. For comparison, optical gain is used as a figureof-merit, and sheet doping density N d 2D or sheet carrier density N s vary from 1.0 × 10 10 cm −2 to 1.0 × 10 12 cm −2 . Doping window (a range of sheet doping/carrier densities with which gain peaks appear) is also confirmed. The calculation results of Pattern A are a reference to other conditions. All calculations are implemented under the heat-  To identify the major causes for the calculation results, the influence of doping conditions on these parameters is investigated. (The plotted dot data are abstracted at the gain peak conditions, and wiggly lines are trend curves.) Furthermore, we also analyze three carrier transition rates determinative for the population inversion: tunneling injection to the upper lasing state, the non-radiative transition between the lasing states, and carrier extraction from the lower-lasing state.
The main component of carrier injection to the upper lasing state is the tunneling transition. When the tunneling between two quantum states ñ |L and ñ |R is considered, the current density of it is described by Eq.

Doping position dependence
The influence of doping positions on optical gain is studied by using narrow-rectangular doping distributions (A-D). As described in Fig. 6, optical gain peaks and doping windows are clearly dependent on doping positions. Among these four patterns, when phonon-wells are doped (A), optical gain takes the highest maximum 63.1 cm −1 , and gain peaks are observed up to around 1.0 × 10 12 cm −2 . However, when doping positions get closer to collection-wells, device performance is exacerbated remarkably. Optical gain of radiation-barrier doping (B) takes almost the same values as the ones for phonon-well doping (A) in low doping conditions. These calculation results, to some extent, support an experimental result implemented by Wang et al. 58) despite small structural differences between two-well and three-well devices. However, the highest gain peak of this doping pattern is 43.0 cm −1 at 4.5 × 10 10 cm −2 , and its doping window shrinks to 2.0 × 10 11 cm −2 . As for other conditions (C and D), the device performance is inferior to these patterns. Comparison in current density and lasing frequency at gain peaks are described in Figs. 7 and 8, respectively. With sheet doping density increasing, current density rises linearly, and  To analyze this result, componential parameters of optical gain; oscillator strength O , 32 population inversion r D , 32 and optical linewidth G 32 are carefully investigated (Fig. 9). All plotted data are abstracted at the fields when optical gain peaks (25-28 kV cm −1 ). As seen in Fig. 9(a), oscillator strength increases when doping positions are proximate to collection-wells, and this behavior is effective to improve optical gain. However, the decrease of population inversion described in Fig. 9(b) surpasses the increase of oscillator strength and determines the behavior of optical gain. Therefore, the simulation results described in Fig. 6 are mainly attributed to the degradation of population inversion. In all doping conditions, population inversion for phononwell doping takes the largest values and gradually decreases from 15% at 1.0 × 10 10 cm −2 to below 5% at 7.0 × 10 11 cm −2 with an increase in doping density. In other doping conditions, on the other hand, population inversion more rapidly decreases. For oscillator strength, the difference due to doping positions is thought to originate from bandbending. Even though band-bending is not so significant in lower doping conditions, wavefunction and eigen energies are very sensitive to potential profile deformation, and their change is reflected in oscillator strength. Optical linewidth increases with doping density except for radiation-barrier doping even though the change is not very impactful on optical gain.
As for the behavior of population inversion, a detailed analysis is implemented. In Fig. 10, the calculated tunneling injection rate ¢ U , 1 3 non-radiative scattering rate between lasing subbands W , 32 and carrier extraction rate W 21 are described where non-radiative scattering rates are the most susceptive to doping positions and determinative to    population inversion. The non-radiative scattering rate of phonon-well doping (A) is around 1.0 × 10 12 s −1 over a whole doping density. However, when the doping position is closer to collection-wells, the non-radiative scattering drastically increases due to excessive ionized-impurity scattering. Moreover, a decrease in the carrier injection rate with doping density also determines the decrease of population inversion. This reduction in injection rate is considered attributable to subband detuning by band-bending effect. Carrier extraction rates are almost immune to these doping conditions.

Doping distribution dependence
The influence of doping distributions and an undoped condition is investigated by comparing the calculation results of the δ-doping (E), homogeneous-doping (F), all layer doping (G) and undoping (H) conditions with those of conventional narrow-rectangular doping (A). As seen in the previous section, maxima of optical gain peaks appear within a doping range 1.0 × 10 10 ∼ 1.0 × 10 12 cm −2 in any doping patterns. First, the δ-doping (E), and homogeneous-doping (F) in phonon-wells are compared with the rectangular doping (A). As observed in Figs. 11-13, doping distributions also influence electrical and optical characteristics greatly. In Fig. 11(a), the optical gain of homogeneous-doping takes 76.2 cm −1 at 3.0 × 10 11 cm −2 which is the second-highest among all patterns, and the one of δ-doping takes almost similar values to the rectangular doping. The doping window     The Japan Society of Applied Physics by IOP Publishing Ltd is observed to become wider when the doping region is wider. Furthermore, under the high doping conditions, a clear difference is seen in current density (Fig. 12) and lasing frequency (Fig. 13), as well. Especially, a discrepancy in current density becomes remarkable, and it reaches higher than 10 kA/cm −2 .
To investigate such behavior in optical gain, the componential parameters are analyzed. In Fig. 14, the difference in these componential parameters amongst the three doping distributions (A, E and F) is not so significant but is surely reflected in the calculation results in Fig. 11. When a doping region is wider, higher oscillator strength and narrower optical linewidth are obtained. Thus, optical gain is also larger in the homogeneous-doping than that in the δ-doping. In the population inversion described in Fig. 14(b), a crossover is seen between the δ-doping and homogeneous-doping at around 2.5 × 10 11 cm −2 . This result is explained by the carrier transition rates shown in Fig. 15. In the low doping conditions, non-radiative scattering rate of homogeneousdoping is higher, and the other two transition rates do not have a so large difference, resulting in a higher population fraction of ULS and population inversion in δ-doping. In the high doping conditions, the carrier injection rate and extraction rate of the homogeneous-doping are also higher than the ones of δ-doping, so the ULS population fractions increases, and the one of LLS decreases, resulting in a higher population inversion with homogeneous-doping where the nonradiative transition rate and carrier extraction rate are enhanced by higher ionized-impurity scattering due to wider doping. To analyze these results based on the behavior of wavefunctions, the band-bending effect induced by localized space charges (electrons and dopants) is observed. In Fig. 16, calculated potential profiles of doping patterns A, F and H with 2.0 × 10 11 cm −2 at 25 kV cm −1 are described. The top and bottom figures describe conduction band edges and Hartree potentials, respectively where different potential shapes are seen depending on the doping distributions. In rectangular doping (A), a concave potential appears in the center of a phonon-well. In homogeneous-doping (F), an almost flat/ideal potential is obtained. In undoped condition (H), where the potential is determined only by electrons, a large difference between the calculated potential and the flat one is seen. In two-well resonant-phonon devices, substantial numbers of electrons accumulate in phonon-wells, resulting in a convex potential as seen in Fig. 16(c). On the other hand, ionized-impurities cause a concave potential. Therefore, doping in phonon-wells cancels out the convex potential created by electrons, and the potential distortion in a whole structure is mitigated.
With impurity doping, band-bending varies depending on doping distributions. As described in Fig. 1, the wavefunction and eigen energy calculated by the Schrödinger equation also shift, along with the complicated potential deformation, and ultimately modify electrical and optical characteristics as described in Figs. 11-13 via componential parameters. Because unique doping distributions generate corresponding   Hartree potentials, the wavefunction and eigen energies also vary. Under the homogeneous-doping condition, the error between the ideal potential and modified potential is the smallest. Thus, detuning of the quantum states is minimized, and high injection and extraction rates are maintained, resulting in high optical gain even at high doping conditions. Practically, the band-bending effect becomes remarkable under high doping conditions, and designed alignment fields in simulation and actual ones confirmed in the experiment are thought to make a large difference. Therefore, homogeneousdoping is expected to solve this problem.
Second, all layer doping (G) is not suitable for enhancing optical gain and the highest operation temperatures. Even though it is advantageous that lasing frequency is immune to doping amount, optical gain is also smaller than that of the reference conditions (A) due to population inversion degradation for the same reason as the doping patterns B-D.
Finally, the undoped condition (H) recorded the highest optical gain 222.7 cm −2 at 4.5 × 10 11 cm −2 , as described in Fig. 11(b), because a rather high oscillator strength and low optical linewidth are obtained even under the high injection condition (Fig. 14). On the one hand, it is reasonable that the narrow optical linewidth is attributed to the absence of ionized-impurity scattering. On the other hand, the increase of oscillator strength is originally attributed to band-bending. With band-bending, a wavefunction overlap of ULS and LLS becomes greater, resulting in higher oscillator strength. The change of the oscillator strength in Fig. 9(a) and nonradiative scattering rate in Fig. 10(b) for doping pattern B-D is also explained by the same idea. The population inversion for the undoped condition is lower than that for the other three conditions (A, E and F) because high tunneling injection works to increase ULS's population but the increased non-radiative scattering rate reduces ULS's population, and the lower carrier extraction rate increases LLS's population. Therefore, population inversion of the undoped condition is degraded.

Discussion
From the simulation results in Sect. 3, the undoped condition is found to be the most preferable setup for high optical gain and high operation temperature in two-well resonant-phonon THz QCLs if only one module is focused. This setup is reasonable when we see that doping regions are separated from lasing regions in MIR QCLs. In reality, however, it is infeasible to drive THz QCLs completely without doping. 54) Therefore, homogeneous-doping in phonon-wells is concluded to be the best method when all modules in an active core are doped. Furthermore, because dopants can migrate and intrude into other layers due to heat during epitaxial growth, it is preferable to avoid doping in the vicinities of interfaces between wells and barriers (Fig. 17).
The knowledge obtained in the simulation is beneficial for expanding doping tactics. Conventionally, an active core of THz QCLs is composed of several hundred modules, and all modules are doped. Therefore, when an active core includes a great number of modules, the total impurity per active core increases simultaneously, resulting in high free-carrier absorption and threshold gain. Having said that, in principle, carriers of QCLs are recycled after radiation, so the authors think that not all modules need doping. A modulation doping scheme, with an active core that includes both doped and undoped modules, should be effective for simultaneously achieving high optical gain and low threshold gain (Fig. 18).
In the recent design trend of THz QCLs, devices composed of two wells have often been studied as a way to improve quantum efficiency by using fewer subbands. On the one hand, decreasing the number of quantum wells per module shortens the length of one module and is effective for high optical gain as described in Eq. (20). On the other hand, it does not change that using an approximately 10 μm thick waveguide is preferable for low waveguide loss. Therefore, down-sizing of devices is expected to increase the net doping density in waveguides if all modules are equally doped.  For example, a long device, designed by Williams et al. with a 539 Å of four-well resonant-phonon structure, has 186 modules per 10 μm active core. 72) On the other hand, a short device, designed by Khalatpour et al. with a 269.3 Å of twowell resonant-phonon structure, has 371 modules in the same size of active core as the long device. 12) Therefore, if all modules are doped in the two devices with the same condition such as around 4.5 × 10 10 cm −2 of sheet doping density over 30 Å-wide near the center of phonon-wells, the net doping density in the short device is almost double that of the long device. Therefore, the highest operation temperature of the short device is thought to be suppressed by high threshold gain due to high free-carrier loss instead of high optical gain. One of the solutions for this issue is a modulation doping scheme, and if every two modules are doped, the net doping density can remain constant as described in Fig. 19. Furthermore, undoped modules are realized in the active core, so a higher optical gain and improved temperature characteristics are also expected. Although instability of the system due to the space charge effect by injection is also apprehensive, 73,74) this issue would     The Japan Society of Applied Physics by IOP Publishing Ltd not be so serious unless the doping period is very long because Straub et al. already demonstrated an operational MIR-QCL device by adopting a modulation doping scheme with longer than 100 nm of doping period in the past. 40) Having said that, they did not discuss the influence of modulation doping on the temperature characteristics of QCLs. In THz QCLs, therefore, the validity of modulation doping schemes need to be confirmed experimentally.
Finally, the duration of carrier recycling is also one of our great interests. The device designed by Straub et al. was already demonstrated to operate even if the doping period is longer than 100 nm. 40) On the other hand, devices completely without doping have also been observed not to lase. 54) Thus, we predict that carrier recycling will have a limitation toward doping period, distance, or space charge effect, and the optimal modulation doping period needs to be adjusted by experiment.

Conclusion
Through theoretical investigation, doping profiles (amount, position, and distribution) have been found to greatly influence the performance of a two-well resonant-phonon THz QCL. In a comparison with a doped single module, we confirmed that impurities should be doped as far as possible from collection-wells where electrons of ULS accumulated. In this case, optical gain is enhanced with doping amount until 2.0-4.0 × 10 10 cm −2 of sheet doping density. Furthermore, wide-doping in phonon-wells cancel out space charge effect (band-bending) due to electrons' accumulation, resulting in a high optical gain. On the other hand, if electrons are able to be injected from an adjacent module to an undoped module, optical gain for the undoped module will be very high due to sharp optical linewidth. This modulation doping scheme, already demonstrated in a MIR-QCL study, is expected to be applied to THz QCLs for enhancement of the highest operation temperature. Even though the detuning of aligned subbands greatly increases under high doping conditions, combining the doping scheme and structural adjustment 63) is expected to further improve device performance.