Simulation of nonequilibrium effects in quantum cascade lasers

Quantum cascade lasers exploit optical transitions between quantized electronic energy levels in multi‐quantum‐well structures for light generation and detection. This design principle enables the realization of compact, semiconductor‐based lasers in the mid‐infrared and terahertz spectral regions. In this contribution, the modeling of such devices is discussed based on an ensemble Monte Carlo approach, enhanced by modifications obtained from the density matrix formalism. A special emphasis is put on the inclusion of nonequilibrium electron and phonon effects. The simulations are validated against available experimental data for a terahertz quantum cascade laser.


Introduction
Quantum cascade lasers (QCLs) are nanostructured semiconductor lasers that exploit intersubband transitions, i.e., optical transitions between quantized energy levels within the conduction band of a multi-quantum-well structure [1]. The lasing wavelength is thus determined by quantum engineering, rather than being dictated by the bandgap energy. In this way, large portions of the terahertz and mid-infrared spectral regions can be covered, ideally complementing the wavelength range that is provided by conventional semiconductor lasers. Also the nonlinear optical properties of the QCL nanostructure are directly related to the quantum design and enable the implementation of novel functionalities. An example is the generation of frequency combs, i.e., broadband comb-like spectra used for sensing and metrology applications.
For a systematic device analysis and design optimization, reliable and computationally efficient simulation approaches are required. In this contribution, stationary carrier transport modeling of QCLs is discussed. Advanced stationary carrier transport simulation approaches, such as the ensemble Monte Carlo (EMC) and nonequilibrium Green's functions (NEGF) methods, treat the relevant scattering mechanisms, for example electron-phonon, electron-impurity and electron-electron interactions, microscopically by evaluating the corresponding Hamiltonians [2]. While NEGF also accounts for coherent transport, EMC is based on stochastic sampling of the Boltzmann transport equation, which constitutes the semiclassical limit of the density matrix (DM) formalism [3]. The restrictions associated with this semiclassical framework can partly be overcome by including corrections derived from the DM approach [3,4]. This allows for example an appropriate treatment of tunneling transport across thick barriers [5], which is indispensable especially for the realistic simulation of various terahertz QCL designs [4]. The resulting DM-EMC approach is also suitable for providing scattering and dephasing rates as input for time dependent DM simulations [6,7], which are used for modeling instability phenomena as well as mode-locked pulse and frequency comb operation of QCLs [8].
Besides its relative numerical efficiency as compared to full quantum transport approaches such as NEGF, a further strength of EMC/DM-EMC is that both non-equilibrium electron and phonon effects can self-consistently be taken into account. This considerably increases the accuracy and versatility of this approach for QCL modeling. In the following, we study the carrier transport in an exemplary experimental terahertz QCL structure, validating our simulations against available experimental results and investigating the influence of nonequilibrium electrons and phonons.

Carrier transport simulation
In the following, we present DM-EMC simulation results for a photon-phonon terahertz QCL design with two wells per period [9]. In Fig. 1(a), the conduction band profile and relevant energy eigenstates are shown, as obtained by a Schrödinger-Poisson solver [10]. Here, the wave functions of levels 3 and 2, corresponding to the upper and lower laser levels, are represented by thick lines. Depopulation of the lower laser level occurs primarily by scattering due to longitudinal optical (LO) phonon emission to level 1 below, which also serves as injector for the upper laser level of the next period. Here, a tight-binding basis is used, where the wave functions are assumed to be localized within their corresponding periods [4].
The carrier transport is modeled based on DM-EMC, which considers all the relevant scattering mechanisms in the structure, including impurity and interface roughness scattering, interactions with LO and acoustic phonons, as well as electronelectron scattering [2]. Since the energy is only quantized in the growth direction of the structure, the electrons can move freely in the in-plane direction of the quantum wells. Thus, the quantized energy levels in QCLs are also referred to as subbands. The ensemble Monte Carlo method takes into account this additional degree of freedom when evaluating the scattering processes. Thus, rather than relying on a priori assumptions, the simulation self-consistently yields for each level the kinetic electron distribution, which can have an effective temperature significantly exceeding the lattice temperature, or even deviate from a thermal distribution [11][12][13]. Also the presence of nonequilibrium LO phonons in excess of thermal occupation, caused mainly by the LO phonon-induced depopulation of the lower laser level, is considered [13][14][15]. In order to account for the photon-related contributions to carrier transport, the simulation is coupled to an evolution equation for the optical resonator field [16,17]. Furthermore, tunneling through the injection barrier is implemented as an additional pseudo-scattering mechanism [4]. For a self-consistent treatment, the broadening of the tunneling and optical transitions must be modeled, which requires the calculation of pure dephasing rates caused by intrasubband scattering, in addition to the lifetime broadening contribution due to intersubband transitions. Similarly as scattering, pure dephasing is here modeled microscopically based on the corresponding Hamiltonians [18,19].

Discussion of nonequilibrium effects
In Fig. 1(b), the measured current-voltage characteristics [9] is compared to simulations with and without including nonequilibrium phonons, yielding good agreement especially for the latter case. Here, pulsed operation at low duty cycle is considered, i.e., the lattice temperature is expected to be close to the sink temperature of T = 10 K. In Fig. 2(a), the kinetic electron distributions in the lasing regime are shown. For nonequilibrium phonons included, a Maxwellian fit yields electron temperatures of 134 K, 37 K and 122 K for the upper, lower and injection levels. This is consistent with available experimental data for terhertz QCLs of the resonant LO phonon depopulation type, yielding electron temperatures of typically 100 K above T for the upper laser level, and close to T for the lower laser level [12]. Furthermore, highly non-Maxwellian features are visible in particular for the kinetic electron distribution of the lower laser level, since its energetic separation from the ground level is below the LO phonon energy [4]. The equilibrium phonon occupation number is given by where k B denotes the Boltzmann constant and E LO is the LO phonon energy. For GaAs with E LO ≈ 36 meV, we obtain N LO = 7.2 × 10 −19 at T = 10 K. Strong LO phonon emission, e.g., in conjunction with the lower laser level depopulation, leads to a significant generation of nonequilibrium (or hot) phonons, subsequently decaying to acoustic phonons with a lifetime of ∼ 5 ps. In our simulation approach, this is captured by a rate equation model, taking into account phonon emission, absorption and decay. The different phonon localizations and wavevector components in growth direction, associated with the transitions between different subbands, are considered by tabulating and tracking the corresponding phonon occupation numbers separately [14]. In Fig. 2(b), the obtained phonon distributions for the most relevant transitions are shown in the lasing regime. The obtained phonon occupations of up to N LO = 0.11 are consistent with experimental data for terahertz QCLs in the lasing regime [15] and exceed the equilibrium occupation of Eq. (1) by 17 orders of magnitude. As can be seen from Fig. 1(b), the simulated current in the lasing regime is somewhat lower, and also closer to experiment, if hot phonons are taken into account, which we attribute to the increased LO phonon absorption and reduced photon-driven current. This reduction of current density at operating bias due to hot phonons has also been observed in EMC simulations of mid-infrared QCLs [13]. While LO phonon emission is not dramatically increased since it scales with N LO + 1, scattering contributions due to LO phonon absorption, which are directly proportional to N LO , can play a relevant role in the presence of hot phonons even at low lattice temperatures. This counteracts the contribution of LO phonon emission to carrier transport, in particular giving rise to thermal backfilling of the lower laser level, and thus reducing the efficiency of LO phonon depopulation. As a consequence, the optical power in the resonator gets reduced, resulting in a smaller contribution of photon-induced processes to carrier transport. As can be seen from Fig. 2(a), also the kinetic electron energy increases slightly due to hot phonons, resulting in somewhat raised effective electron temperatures (10..20 K) as compared to simulations without hot phonons. For completeness, the contribution of phonons to pure dephasing is also considered in our simulations, enhancing our modeling framework. However, in the example chosen here, pure dephasing is dominated by elastic scattering due to interface roughness and impurities, as has also been found in other theoretical studies of intersubband quantum well structures [19,20]. As mentioned in Section 2, the optical cavity field is self-consistently included in our simulation. However, a direct comparison between simulated and measured optical power is difficult since the collection efficiency of the measurement setup is hard to determine exactly in the terahertz regime. An indirect measure of the optical power is the internal quantum efficiency, extracted by comparing the currents at operating bias of a lasing structure and a non-lasing one, realized in experiment by adding absorption to the waveguide. The difference in currents is then attributed to the photon-driven contribution, and the internal quantum efficiency is given by the ratio of photon-driven to total current [9]. Our simulation yields a lower efficiency (46.5%) when hot phonons are included, which is also closer to the experimental value (43%) [9], than without hot phonons (51.9%). This again confirms the importance of including hot phonon effects for accurate QCL simulations.

Conclusion
In conclusion, we have presented carrier transport simulations of a terahertz QCL structure. Here we have employed a DM-EMC approach which is based on the semiclassical EMC method, but includes a microscopic model of pure dephasing and a treatment of tunneling derived from the density matrix formalism. A main advantage of EMC is that it calculates the kinetic electron distributions in the subbands, and thus implicitly considers nonequilibrium electron effects. Furthermore, nonequilibrium phonons can be taken into account. Our approach also self-consistently includes the coupling between the carrier transport and the optical cavity field. By comparison of simulation results to available experimental data, we have shown that the proper inclusion of nonequilibrium electron and phonon effects is important to obtain good agreement, hence underlining the importance of such effects in QCLs. In particular, the presence of hot phonons increases LO phonon absorption, resulting in thermal backfilling of the lower laser level. As a consequence, the simulation yields a reduced optical power and a smaller current in the lasing regime, which agree well with experimental data.