Quantum cascade lasers are unipolar semiconductor lasers covering a wide range of the infrared and terahertz spectrum. Lasing action is achieved by using optical intersubband transitions between quantized states in specifically designed multiple-quantum-well heterostructures. A systematic improvement of quantum cascade lasers with respect to operating temperature, efficiency, and spectral range requires detailed modeling of the underlying physical processes in these structures. Moreover, the quantum cascade laser constitutes a versatile model device for the development and improvement of simulation techniques in nano- and optoelectronics. This review provides a comprehensive survey and discussion of the modeling techniques used for the simulation of quantum cascade lasers. The main focus is on the modeling of carrier transport in the nanostructured gain medium, while the simulation of the optical cavity is covered at a more basic level. Specifically, the transfer matrix and finite difference methods for solving the one-dimensional Schrödinger equation and Schrödinger-Poisson system are discussed, providing the quantized states in the multiple-quantum-well active region. The modeling of the optical cavity is covered with a focus on basic waveguide resonator structures. Furthermore, various carrier transport simulation methods are discussed, ranging from basic empirical approaches to advanced self-consistent techniques. The methods include empirical rate equation and related Maxwell-Bloch equation approaches, self-consistent rate equation and ensemble Monte Carlo methods, as well as quantum transport approaches, in particular the density matrix and non-equilibrium Green's function formalism. The derived scattering rates and self-energies are generally valid for n-type devices based on one-dimensional quantum confinement, such as quantum well structures.

## I. INTRODUCTION

The quantum cascade laser (QCL) is a fascinating device, combining aspects from different fields such as nanoelectronics and quantum engineering, plasmonics, as well as subwavelength and nonlinear photonics. Since its first experimental realization in 1994,^{1} rapid progress has been achieved with respect to its performance, in particular the output power and the frequency range covered. The development of innovative types of QCLs and subsequent design optimization has gone hand in hand with detailed modeling, involving more and more sophisticated simulation tools. The advancement of simulation methods is driven by an intrinsic motivation to describe the underlying physical processes as exactly and realistically as possible. In addition to scientific motivations, improved agreement with experimental data and enhanced predictive power is desired. Furthermore, the extension and improvement of QCL simulation methods is largely driven by experimental progress, such as innovative designs.

The QCL is a special type of semiconductor laser where the optical gain medium consists of a multiple quantum well heterostructure, giving rise to the formation of quantized electron states. These so-called subbands assume the role of the laser levels in the QCL. To date, lasing has only been obtained for n-type QCLs, i.e., designs using intersubband transitions in the conduction band. Well established material systems are InGaAs/InAlAs on InP substrate, and GaAs/AlGaAs on GaAs.^{2,3} Furthermore, antimonide QCLs have been demonstrated, e.g., near- and mid-infrared QCLs using InAs/AlSb^{4,5} and aluminum-free InGaAs/GaAsSb QCLs operating in the mid-infrared^{6} and terahertz^{7} regime. Since all these materials (apart from AlSb) have direct bandgaps, the development of simulation methods and tools for QCLs has focused on the conduction band Γ valley where the lasing transitions take place.

In contrast to conventional semiconductor lasers where lasing occurs due to electron-hole recombination between conduction and valence band, the emission wavelength of QCLs is not determined by the bandgap of the material. Rather, by adequately designing the nanostructured gain medium, the device can be engineered to lase at a certain frequency. In this way, the mid/far infrared and terahertz regimes become also accessible, important for applications such as chemical and biological sensing, imaging, and special communication applications. In addition, the QCL offers the typical advantages of semiconductor devices, such as compactness, reliability, and potentially low production costs.^{8} The basic working principle of the QCL is illustrated in Fig. 1. Lasing is obtained due to stimulated emission between the upper and the lower laser level, here resulting in mid-infrared radiation at a wavelength of 5 *μ*m. The depopulation level is separated from the lower laser level by the longitudinal optical (LO) phonon energy for InGaAs of ∼30 meV, enabling efficient depopulation of the lower laser level. The injector level, located between two adjacent periods, ensures electron transport and injection into the upper laser level of the next period. The optical field confinement and beam shaping is provided by a specifically designed resonator, which is frequently based on plasmonic effects and can exhibit subwavelength structuring for a further performance improvement.

The QCL operating principle was devised in 1971 by Kazarinov and Suris.^{10} However, due to the experimental challenges involved, it took more than two decades until QCL operation was experimentally demonstrated at Bell Labs by Capasso and his co-workers and Faist *et al.*^{1} Since then, considerable progress has been achieved, e.g., with respect to operating temperature,^{11} output power,^{12} and the available frequency range.^{5,13} Meanwhile, QCLs cover a spectral range from 2.6 *μ*m to above 400 *μ*m (obtained by applying an external magnetic field).^{5,13} At room temperature, continuous wave output powers of several W can now routinely be obtained in the mid-infrared region,^{12} and widely tunable designs are available.^{14} Furthermore, at cryogenic temperatures wall-plug efficiencies of 50% and above have been demonstrated,^{9,15} and up to 27% at room temperature.^{12} Commercial QCLs and QCL-based systems are already available from various companies.^{8}

On the other hand, terahertz QCLs, first realized in 2002,^{16} still require cryogenic operating temperatures. Lasing up to ∼200 K has been obtained,^{17} and by applying a strong external magnetic field, an operating temperature of 225 K has been demonstrated.^{18} A long standing goal is to achieve operation at room temperature or at least in the commercially available thermoelectric cooler range (∼240 K). For a further enhancement of the operating temperature, it is crucial to improve the gain medium design and reduce the resonator loss.^{17,19} For systematic QCL optimization and exploration of innovative designs, detailed modeling is essential, requiring advanced simulation approaches which enable realistic QCL simulations.

Besides its function as the optical gain medium, the QCL quantum well heterostructure can also exhibit a strong nonlinear optical response. Like the laser transition itself, such nonlinearities are based on quantized electron states, and thus can be custom-tailored for specific applications and optimized to exhibit extremely high nonlinear susceptibilities. This has led to a new paradigm for QCL-based THz generation at room temperature: Here, difference frequency mixing is used, where the THz radiation corresponds to the difference in frequency of two detuned mid-infrared beams, generated by highly efficient mid-infrared QCLs.^{20} The optical nonlinearity is either implemented into the QCL separately in a specifically designed heterostructure,^{21} or, more commonly, directly integrated into the gain medium.^{20} With such an approach, THz generation at room temperature was obtained at power levels of up to 120 *μ*W.^{22,23} Furthermore, broadband frequency tunability has been demonstrated, which is important for various applications such as spectroscopy.^{23,24} Here, the main goal is to push the available room temperature output power to a few mW, as required for most technical applications. The modeling of such terahertz QCL sources is particularly demanding, requiring the coupled simulation of two mid-infrared QCLs and a careful modeling of the nonlinear susceptibility and the associated frequency conversion process in the heterostructure.^{25} Artificial optical nonlinearities have also been used to extend QCL operation towards shorter wavelengths, which cannot be reached directly since the energy spacing of the quantized laser levels is limited by the quantum well depth of the gain medium. Based on a frequency doubling structure,^{26} a room temperature QCL source operating at 2.7 *μ*m could be demonstrated.^{27}

The ongoing development and optimization of the QCL heterostructure is accompanied by improvements of the optical cavity. For example, changes in the design and material of the plasmonic resonator structure have played a crucial role for increasing the operating temperature of THz QCLs.^{17,19} Also, the performance of difference frequency THz sources has largely benefited from special cavity designs, such as the surface emitting Cherenkov waveguide scheme to obtain increased output power and efficiency.^{23,28} Thus, a quantitative modeling and reliable numerical design optimization and exploration must also include the optical cavity. Moreover, various types of resonators with periodic subwavelength structuring have been developed in the THz and mid-infrared regime. This includes distributed feedback structures to enforce single mode lasing and distributed Bragg reflectors for enhanced facet reflectivity.^{29–31} A further example are surface emission schemes based on one- and two-dimensional photonic crystal structures,^{29,32,33} offering tailorable emission properties and improved beam quality. The simulation of such subwavelength-structured cavities requires advanced electromagnetic modeling, e.g., based on coupled mode theory^{34} or even full finite difference time domain simulations of Maxwell's equations.^{33,35}

The goal of this review is to give a detailed survey and discussion of the modeling techniques used for QCL simulation, ranging from basic empirical approaches to advanced self-consistent simulation methods. The focus is here on the modeling of the heterostructure gain medium. Also the simulation of the optical cavity will be covered for simple resonator waveguide structures. As mentioned above, the modeling of complex cavities such as photonic crystal structures constitutes an advanced electromagnetic modeling task, which is beyond the scope of this paper. The review is organized as follows: In Sec. II, the numerical solution of the one-dimensional Schrödinger equation is discussed, providing the eigenenergies and wave functions of the energy states in the QCL heterostructure. Furthermore, the inclusion of space charge effects by solving the Schrödinger-Poisson equation system is treated. Section III covers the modeling of the optical resonator, where we focus on a basic waveguide resonator structure. In Sec. IV, an overview and classification of the different carrier transport models is given which are commonly used for the theoretical description of the QCL gain medium. Section V contains a discussion of empirical modeling approaches for the gain medium, relying on experimental or empirical input parameters, namely, empirical rate equation and related Maxwell-Bloch equation approaches, which are discussed in Secs. V A and V B, respectively. Advanced self-consistent methods, which only require well known material parameters as input, are covered in Secs. VI–VIII: Semiclassical approaches such as the self-consistent rate equation model (Sec. VI E) and the ensemble Monte Carlo method (Sec. VI F), as well as quantum transport approaches such as the density matrix method (Sec. VII) and the non-equilibrium Green's function formalism (Sec. VIII). In this context, also the transition rates (Secs. VI A–VI D) and self-energies (Sec. VIII D) are derived for the relevant scattering processes in QCLs. The paper is concluded in Sec. IX with an outlook on future trends and challenges.

## II. SCHRÖDINGER-POISSON SOLVER

A careful design of the quantized states in the QCL heterostructure is crucial for the development and optimization of experimental QCL devices. In particular, the lasing frequency is determined by the energy difference between the upper and lower laser level. Furthermore, a careful energy alignment of the levels is necessary to obtain an efficient injection into the upper and depopulation of the lower laser level. For example, for the structure shown in Fig. 1, the depopulation level is separated from the lower laser level by the InGaAs LO phonon energy of ∼30 meV to obtain efficient depopulation, while the injector level is aligned with the upper laser level of the next period.^{9} Also careful engineering of the wave functions is important, determining the strength of both the optical and nonradiative transitions. Moroever, various carrier transport simulation approaches, such as ensemble Monte Carlo (EMC), require the quantized energy levels as input. The eigenenergies and wave functions are determined by solving the stationary Schrödinger equation or, if space charge effects are taken into account, the Schrödinger-Poisson equation system.

### A. Schrödinger equation

The QCL heterostructure consists of alternately grown thin layers of different semiconductor materials, resulting in the formation of quantum wells and barriers in the conduction band. Thus, the potential *V* and effective masses vary in the growth direction, here denoted by *z*, whereas *x* and *y* refer to the in-plane directions. The QCL heterostructure is usually treated in the framework of the Ben Daniel-Duke model which only considers the conduction band,^{36} where the lasing transitions and carrier transport take place. Within this approximation, the stationary Schrödinger equation is given by

where $\psi 3D$ and *E* denote the wave function and eigenenergy, respectively. Furthermore, $m\u2225$ refers to the in-plane effective mass and *m** is the effective mass in the growth direction, i.e., perpendicular to epitaxial layers. For bound states, the wave function is commonly normalized, i.e., $\u222c\u222b\u2212\u221e\u221e|\psi 3D|2dxdydz=1$. Since *V* and the effective masses only depend on the *z* coordinate, we can make the ansatz

Here, *S* is the in-plane cross section area and $k=[kx,ky]T$ denotes the in-plane wave vector, where T indicates the transpose. The factor *S*^{−1∕2} is added in Eq. (2) to obtain the normalization condition $\u222b|\psi k|2dz=1$. Insertion of Eq. (2) in Eq. (1) yields the Ben Daniel-Duke model

where the wave function $\psi k(z)$ and energy *E*** _{k}** depend on the in-plane electron motion, i.e., on

**k**. Decoupling can be obtained by neglecting the

*z*dependence of the in-plane effective mass $m\u2225$, yielding the one-dimensional (1D) Schrödinger equation in its usual form

The eigenenergy *E* due to 1D electron confinement in *z* direction is related to the total energy *E*** _{k}** by

*E*

**=**

_{k}*E*+

*E*

_{kin}, where

is the kinetic electron energy due to the free in-plane motion of the electrons. The effective mass values not only depend on the material composition^{37} but also on the lattice temperature and doping level.^{38} The latter effects tend to play a secondary role in QCLs and are thus usually neglected. For strained QCL structures, the effective masses are additionally affected by the lattice mismatch between the different semiconductor materials, resulting in modified values $m\u2225\u2260m\u2217$.^{39}

### B. Boundary conditions

Strictly speaking, the quantum states in the QCL heterostructure are not bound, since the electron energy *E* exceeds the barrier potential for large values of *z*. Thus, the electrons will tunnel out of the multiple quantum well system after a limited time. This situation is illustrated in Fig. 2. If the electron remains in the quantum well system for a considerable amount of time, the concept of quasi-bound or quasi-stationary states can be used.^{40} Such a state can again be described by Eq. (4), however the eigenenergy is now a complex quantity, $E=E0\u2212i\u210f\gamma q/2$. The meaning of the imaginary term $\u2212i\u210f\gamma q/2$ for quasi-stationary states can be understood by considering that the time dependence of a stationary wave function ψ is given by $exp(\u2212iEt/\u210f)$. Consequently, for a complex eigenenergy *E* the density probability $|\psi |2$ of the electron in the quantum well system decays as $exp(\u2212\gamma qt)$, thus $\gamma q$ is the probability density decay rate of the quasi-bound state, and $\tau q=\gamma q\u22121$ can be associated with the lifetime of the particle in the quantum well.^{40}

Some numerical approaches have been developed to solve Eq. (4) for quasi-bound states.^{41,42} However, the relevant states which are considered for the design of a QCL structure, such as the upper and lower laser levels and injector states, are typically strongly bound to obtain optimum performance. Since the electron lifetime in these states is governed by scattering processes and resonant tunneling rather than $\tau q$, they are usually treated as bound states. This can be done by restricting the simulation to a finite simulation window containing a limited number of periods, and imposing artificial boundary conditions $\psi =0$ at the borders, as illustrated in Fig. 3. The simulation window has to be chosen sufficiently large so that the portion of the wave function which lies inside the quantum wells, i.e., which significantly deviates from 0, is not markedly affected. For the simulation, we assume that the QCL heterostructure is periodic, where a period is defined by a sequence of multiple barriers and wells, as illustrated in Fig. 1. Thus, for a biased QCL structure, it is sufficient to compute the eigenenergies and corresponding wave functions for a single energy period of width *E*_{p}, corresponding to the bias drop over a QCL period (see Fig. 3). The solutions of the other periods can then simply be obtained by shifts in position by multiples of the period length *L*_{p}, and corresponding shifts in energy. For example, the solutions in the right-neighboring period are given by $\psi i\u2032(z)=\psi i(z\u2212Lp),\u2009Ei\u2032=Ei\u2212Ep$, where $\psi i$ and *E _{i}* are the wave function and eigenenergy of the

*i*th solution of Eq. (4) obtained in the central period.

For the laser design and simulation, usually only the relevant laser subbands such as the upper and lower laser level and injector states are considered. These typically correspond to the strongly bound states of the QCL heterostructure. Furthermore, the artificial boundary conditions illustrated in Fig. 3 give rise to spurious solutions of Eq. (4), which are not localized in the quantum wells. A systematic selection of the relevant subbands can be achieved by considering only the most strongly bound levels in each period. In this context, the energy of state *i* relative to the conduction band edge

is a useful quantity.^{43} An automated selection of the relevant subbands can then be implemented by choosing the states with the lowest energies $E\u0303i$ in each period, considering all subbands which contribute significantly to the carrier transport.

### C. Nonparabolicity

In Fig. 4, the band structure of GaAs as a typical III–V semiconductor material is displayed. Shown is the valence band (dashed lines), consisting of heavy hole, light hole, and split-off bands, and the conduction band (solid line). As pointed out above, only the conduction band is considered in the Ben Daniel-Duke model, introduced in Sec. II A. The conduction band has three minima, referred to as the Γ, L, and X (or Δ) valley according to their position in **K** space. For direct bandgap semiconductors typically used in QCLs, the Γ valley is the lowest minimum. Thus, the theoretical treatment is usually restricted to that valley, although under certain conditions transitions to the X or L valleys can also affect QCL operation.^{45–48} Furthermore, Eqs. (1)–(5) assume a parabolic dispersion relation between energy and wave number (dotted line in Fig. 4). This approximation only holds close to the Γ valley minimum, i.e., it breaks down for high-lying energy levels. The deviation from the parabolic dispersion, i.e., the nonparabolicity, scales roughly inversely with the semiconductor bandgap.^{49} This effect thus plays a role especially for mid-infrared QCLs with lasing transition energies of up to a few 100 meV, which are frequently based on low-bandgap semiconductors such as InGaAs or InGaSb. A detailed treatment of nonparabolicity effects in QCLs has been performed by using **k·p** theory.^{45–47,50} A related strategy, which we will discuss in the following, is to consider nonparabolicity in Eq. (4) by implementing an energy dependent mass, derived from **k·p** theory. This method involves additional approximations, but is less complex than full **k·p** calculations and can easily be incorporated into the transfer matrix approach, which is widely used to solve Eq. (4). According to Ekenberg,^{51} we obtain the energy dependent effective masses

Here, $\alpha \u2032$ and $\beta \u2032$ are nonparabolicity parameters, which are defined in the framework of a 14-band **k·p** calculation.^{51,52} Approximately, we have $\alpha \u2032=(Eg+\Delta so/3)\u22121$,^{49} where *E*_{g} and $\Delta so$ are the bandgap energy and the energy difference between the light-hole and split-of valence bands, respectively. Furthermore, *m** is the Γ valley effective mass. In the QCL heterostructure, the electron kinetic energy in the conduction band is given by $E\u2032(z)=E\u2212V(z)$. Since the semiconductor material changes along the growth direction, also $m\u2217,\u2009\alpha \u2032$, and $\beta \u2032$ are *z* dependent. This approach is valid for moderate energies $E\u2032\u226a1/\alpha \u2032$. The fact that Eq. (7a) even breaks down for $\alpha \u2032E\u2032>1/4$ can cause numerical problems for large values of *z* close to the right simulation boundary, where ψ is already close to 0 but $E\u2032$ can assume relatively large values (see Fig. 3). This issue can be avoided by using Eq. (7a) for $E\u2032<0$, and its second order Taylor expansion

for $E\u2032>0$, which is the widely used lowest order implementation of nonparabolicity.^{49,51} For an energy dependent effective mass, the Hamiltonian in Eq. (4) is not Hermitian; thus, the obtained wave functions are in general not orthogonal.^{53} Also other approaches are available for implementing nonparabolicity effects, for example based on the Kane model.^{53}

In self-consistent simulations, each scattering process is evaluated based on the corresponding Hamiltonian, which depends on the in-plane effective masses of the subbands involved. Thus, if we want to include corrections due to nonparabolicity also for the scattering processes,^{54,55} then Eq. (7b) should be used to determine $m\u2225(E\u2032)$, which is however *z* dependent. This problem can be avoided by defining an average effective mass for the *i*th subband

where the wave function is assumed to be normalized, $\u222b|\psi i|2dz=1$, and *E _{i}* denotes the corresponding subband eigenenergy. This lowest order implementation of nonparabolicity thus yields different, albeit constant in-plane masses for each subband.

### D. Numerical solution

Numerical approaches for solving the one-dimensional effective mass Schrödinger equation [Eq. (4)] are required to be robust. Also computational efficiency is crucial, especially for QCL design and optimization tasks where many simulations have to be performed. Furthermore, a straightforward implementation is desirable. Widely used numerical approaches include the transfer matrix method^{41,43,56,57} and finite difference scheme.^{58,59} Both methods have their strengths and shortcomings. In particular, effects such as nonparabolicity can be included more easily into the transfer matrix approach. A further advantage is the exact treatment of the potential steps between barriers and wells in the QCL heterostructure. On the other hand, this method can exhibit numerical instabilities for multiple or extended barriers due to an exponential blowup caused by roundoff errors.^{58} This issue can however be overcome, for example by using a somewhat modified approach, the scattering matrix method.^{60} In Fig. 5, the transfer matrix and finite difference schemes are illustrated.

#### 1. Transfer matrix approach

The transfer matrix approach uses the fact that analytical solutions of Eq. (4) are available for constant or linear potential sections and for potential steps.^{61} An arbitrary potential can then be treated by approximating it in terms of piecewise constant or linear segments, respectively. In the first case, the solution is given by complex exponentials,^{41,56} while the approximation by linear segments gives rise to Airy function solutions.^{41} The Airy function approach provides an exact solution for structures which consist of segments with constant effective masses and piecewise linear potentials, such as a biased QCL heterostructure if we neglect nonparabolicity and space charge effects. On the other hand, Airy functions are much more computationally expensive than exponentials, and also prone to numerical overflow for segments with nearly flat potential.^{62} Thus, great care has to be taken to avoid numerical problems and to evaluate the Airy functions efficiently.^{63}

In the following, we focus on the exponential transfer matrix scheme, illustrated in Fig. 5(a). We start by dividing the structure into segments which can vary in length, but should be chosen so that all band edge discontinuities occurring between wells and barriers are located at the border of a segment, i.e., do not lie within a segment. The potential and effective mass in each segment *n* with $zn\u2264z<zn+\Delta n=zn+1$ are approximated by constant values, e.g., $Vn=V(zn),\u2009mn\u2217=m\u2217(zn)$, resulting in a jump $Vn\u2192Vn+1,mn\u2217\u2192mn+1\u2217$ at the border between the segments *n* and *n* + 1.^{41} Nonparabolicity can straightforwardly be implemented by using Eq. (7a) for *m**; then $mn\u2217=mn\u2217(E)$ depends on the eigenenergy *E*. The solution of Eq. (4) in segment *n* is given by

Here, $kn=2mn\u2217(E\u2212Vn)/\u210f$ is the wave number (for $E<Vn$, we obtain $kn=i\kappa n=i2mn\u2217(Vn\u2212E)/\u210f$).^{41} The matching conditions for the wave function at a potential or effective mass discontinuity read

where *z*_{0}+ and *z*_{0}– denote the positions directly to the right and left of the discontinuity.^{61} Using Eqs. (10) and (11), the amplitudes *A _{n}*

_{+1}and

*B*

_{n}_{+1}can be related to

*A*and

_{n}*B*by

_{n}where the transfer matrix is with $k\u0303n=kn/mn\u2217$ given by

Here,

corresponds to the transfer matrix for a flat potential obtained from Eq. (10), and

is the potential step matrix for the interface between the segments *n* and *n* + 1 at *z*_{0} = *z _{n}*

_{+1}, derived from Eq. (11).

^{61}For a structure divided into

*N*segments, the relation between the amplitudes at the left and right boundaries of the structure,

*A*

_{0},

*B*

_{0}and

*A*,

_{N}*B*, can be obtained by multiplying the matrices for all segments

_{N}This equation must be complemented by appropriate boundary conditions. As described in Sec. II B, for the QCL heterostructure we can restrict our simulation to a limited number of periods (see Fig. 3) and assume $\psi =0$ at the boundaries of our simulation window, corresponding to *A*_{0} + *B*_{0} = 0 and *A _{N}* +

*B*= 0. The left boundary condition can, for example, be enforced by setting

_{N}*A*

_{0}= 1,

*B*

_{0}= −1. The right boundary condition

*A*+

_{N}*B*= 0 can then only be satisfied if the energy dependent matrix elements

_{N}*T*

_{11}(

*E*),

*T*

_{12}(

*E*),

*T*

_{21}(

*E*), and

*T*

_{22}(

*E*) in Eq. (16) assume certain values. The corresponding energies

*E*are the eigenenergies of the bound states.

^{41}Numerically, these eigenenergies are found by the so-called shooting method. Here, the wave function at the right boundary is computed from Eq. (16) as a function of energy, $\psi N(E)=AN(E)+BN(E)$, and the eigenenergies are given by $\psi N(E)=0$. For the periodic QCL heterostructure, it is sufficient to restrict the simulation to a single energy period (see Fig. 3). In practice, Eq. (16) can only be solved for a limited number of discrete energy points $Em=E0+m\Delta E$ with sufficiently close spacing $\Delta E$. The eigenenergies are located in intervals $Em..Em+1$ with $\psi N(Em)\psi N(Em+1)<0$, and can be determined more accurately by choosing a finer energy grid in the corresponding intervals, or preferably by applying a root-finding algorithm such as the bisection method.

^{64}

The accuracy of the transfer matrix scheme can generally be improved by replacing Eq. (13) with a symmetric transfer matrix

#### 2. Finite difference method

The finite difference method works by converting Eq. (4) to a finite difference equation. As illustrated in Fig. 5(b), a spatial grid with uniform spacing $\Delta z$ is introduced, and the wave function $\psi (z)$, potential $V(z)$, and effective mass $m\u2217(z)$ are represented by the corresponding values $\psi n$, *V _{n}*, and $mn\u2217$ on the grid points

*z*. First order derivatives are approximated by $\u2202z\psi (zn+\Delta z/2)\u2248\Delta \psi n+1/2=(\psi n+1\u2212\psi n)/\Delta z$. Consequently, the term $\u2202z(m\u2217)\u22121\u2202z\psi $ in Eq. (4) at

_{n}*z*=

*z*can be expressed as $(\Delta \psi n+1/2/mn+1/2\u2217\u2212\Delta \psi n\u22121/2/mn\u22121/2\u2217)/\Delta z$. Using linear interpolation $mn+1/2\u2217=(mn\u2217+mn+1\u2217)/2$, we obtain the discretized form of Eq. (4)

_{n}^{58,59}

with

^{58}as illustrated in Fig. 5(b). We again assume that the wave function is zero at the boundaries of our simulation window, $\psi 0=\psi N=0$. Writing Eq. (18) in matrix form yields an eigenvalue equation of the form $(H\u2212EI)\psi =0$, where $\psi $ is the wave function vector with $\psi =[\psi 1,\psi 2,\u2026,\psi N\u22121]T$,

**I**represents the identity matrix of size

*N*– 1, and

**H**is the Hamiltonian matrix with the non-zero elements

*H*

_{n}_{,}

_{n}=

*d*,

_{n}*H*

_{n}_{,}

_{n}

_{–1}= −

*s*,

_{n}*H*

_{n}_{,}

_{n}

_{+1}= −

*s*

_{n}_{+1}. This equation can be solved by using a standard eigenvector solver for a tridiagonal matrix problem.

^{64}However, if nonparabolicity is taken into account, this is not directly possible since the effective masses $mn\u2217$ depend on

*E*, and a modified numerical scheme must be employed.

^{65}Alternatively, the shooting method discussed in Sec. II D 1 can be employed.

### E. Schrödinger-Poisson equation system

To lowest order, electron-electron interaction can be considered by self-consistently solving Eq. (4) together with the Poisson equation^{57,66}

From a quantum mechanical point of view, this corresponds to a mean-field treatment of the electron-electron interaction referred to as Hartree approximation, representing the lowest order of a perturbation expansion in the electron-electron interaction potential. In Eq. (20), $\u03f5(z)$ is the permittivity which varies with semiconductor composition and thus is also periodic, and *e* denotes the elementary charge. The right hand side of Eq. (20) corresponds to the space charge ρ in the QCL heterostructure due to the positively charged donors with concentration $nD(z)$ and the electrons, where $nis$ is the electron sheet density of level *i* with wave function $\psi i(z)$. This charge distribution in the structure gives rise to space charge effects, resulting in an additional electrostatic potential energy $V\u0303(z)$ which causes conduction band bending.^{67}

The total potential *V* in Eq. (4) is then given by $V=V0+V\u0303$. Here, $V0=Vc\u2212Epz/Lp$, where *V*_{c} is the unbiased conduction band profile due to the varying material composition, thus describing the wells and barriers, and the term –*E*_{p} *z/L*_{p} results from the applied bias. Since the energy drop across a period is given by the external bias, $V(z0)\u2212V(z0+Lp)=Ep$, we have $V\u0303(z0)=V\u0303(z0+Lp)$. Due to the charge neutrality in each period, $\u222bz0z0+Lp\rho dz=0$, we furthermore obtain $\u2202zV\u0303(z0)=\u2202zV\u0303(z0+Lp)$, i.e., $V\u0303$ has the periodicity of *V*_{c}. Thus, we can restrict the solution of Eq. (20) to a single QCL period $z\u03f5[z0,z0+Lp]$ and assume the boundary conditions $V\u0303(z0)=V\u0303(z0+Lp)=0$.

Equation (20) can, for example, be solved by applying the finite difference method. In analogy to Eqs. (18) and (19), we obtain

with

and

Equation (21) is then solved over a single QCL period $z\u03f5[z0,z0+Lp]$, with the grid points *z _{n}*, $n=0\u2026P$ which should coincide with the grid used for solving the Schrödinger equation Eq. (4). Applying the boundary conditions $V\u03030=V\u0303P=0$, Eq. (21) can be written as a matrix equation $MV\u0303=\rho $, where $V\u0303$ and $\rho $ represent vectors with elements $V\u0303n$ and $\rho n$, respectively, with $n=1\u2026(P\u22121)$.

**M**is a matrix with the non-zero elements $Mn,n=\u2212d\u0303n,\u2009Mn,n\u22121=s\u0303n,Mn,n+1=s\u0303n+1$. This equation can be efficiently solved using an algorithm for tridiagonal equation systems.

^{64}

While $nis$ in Eq. (22) can in principle only be determined by detailed carrier transport simulations, simpler and much faster approaches are often adopted, e.g., for design optimizations of experimental QCL structures over an extended parameter range. Frequently, Fermi-Dirac statistics is applied,^{57,66} with

Here, we assume that the electron distribution is described by the lattice temperature *T*_{L}. Furthermore, μ is the chemical potential, *k*_{B} denotes the Boltzmann constant, and $mi\u2225$ is the effective mass associated with the *i*th subband, see Eq. (9). If nonparabolicity effects are neglected, $mi\u2217$ can often be approximated by the value of the well material. In Eq. (24), the energy $E\u0303i$ defined in Eq. (6) is used instead of the eigenenergy *E _{i}*, to correctly reflect the invariance properties of the biased structure.

^{43}Especially, this guarantees that the obtained results do not depend on the choice of the elementary period in the heterostructure. The chemical potential μ is found from the charge neutrality condition within one period, i.e.,

This is done recursively by first determining a lower and an upper boundary value for μ where $\u2211inis<ns$ and $\u2211inis>ns$, respectively, and then finding the exact μ, e.g., by using the bisection method.^{64}

The total potential in the Schrödinger equation (4) is given by $V=V0+V\u0303$, where $V\u0303$ has to be obtained by solving the Poisson equation (20). On the other hand, the wave functions $\psi i$ in Eq. (20) must be determined from Eq. (4). In practice, this is done by iteratively solving the Schrödinger and Poisson equations,^{57} initially assuming $V\u0303=0$, until the results for $V\u0303,\u2009\psi i$, and *E _{i}* converge. If the $nis$ are obtained from self-consistent carrier transport simulations instead of using Eq. (24), then the carrier transport simulation and the numerical solution of the Schrödinger-Poisson system have to be performed iteratively until convergence is obtained. Such simulations are then referred to as self-self-consistent approaches.

^{67,68}

In Fig. 6, simulation results for a terahertz QCL^{69} obtained by solving the Schrödinger equation (solid lines) and the Schrödinger-Poisson system (dashed lines) are compared. Although the populations obtained by Eq. (24) usually deviate considerably from detailed carrier transport simulations, this simplified approach can already give improved results, as compared to completely neglecting space charge effects.^{68} The reason is that the positively charged donors are usually localized in a relatively small section, e.g., a single quantum well, while the electrons are more or less distributed across the whole QCL period. This leads to considerable conduction band bending due to space charge effects, as illustrated in Fig. 6. While a more exact determination of $nis$ yields a somewhat different electron distribution, the overall result for $V\u0303$ will be similar as for the simplified model based on Eq. (24).

The method shown above is not the only one for solving the Schrödinger-Poisson system self-consistently under the condition of global charge neutrality. A common alternative approach is to keep the chemical potential fixed and solve the Poisson equation with Neumann boundary conditions. These boundary conditions allow the Poisson potential to self-adjust the density and maintain the global charge neutrality. Typically, this requires a nonlinear realization of the Poisson equation to include an explicit potential dependence of the charge density. Such a nonlinear Poisson equation can then effectively be solved with iterative methods such as the predictor-corrector approach.^{70}

## III. OPTICAL RESONATOR MODELING

The resonator spatially confines the radiation field, furthermore providing optical outcoupling, beam shaping and frequency selection. While most theoretical work has focused on the carrier transport in the gain medium, there has also been progress in the modeling of the cavity. For example, the resonator loss, which is crucial for the temperature performance of THz QCLs,^{17,19} has been extracted from finite element simulations of the resonator.^{71} In Fig. 7, a typical waveguide resonator geometry of a QCL is sketched, and the used coordinate system is shown for reference. The propagation direction of the optical field is denoted by *x*, *y* refers to the lateral direction, and *z* indicates the growth direction of the heterostructure. The resonator consists of materials with different permittivities to obtain waveguiding and optical outcoupling. Thus, for optical cavity simulations, the material permittivities must be known, which in general depend on the frequency, doping level, and temperature. In this context, often the Drude model is employed with adequately chosen fitting parameters.^{71} For intersubband optical transitions, only the dipole matrix element in the *z* direction where quantum confinement occurs is nonzero, see Sec. V B. Thus, only resonator modes with an electric field component along the *z* direction are amplified. For this reason, surface emission, i.e., outcoupling through the *xy*-plane, can only be obtained by using special outcoupling schemes.^{33,72}

### A. Maxwell's equations

The propagation of the electric and magnetic field vectors $E(x,t)$ and $H(x,t)$ is generally described by Maxwell's equations. Here, $x=[x,y,z]T$ and *t* are the position vector and time, respectively. In the following, monochromatic fields are considered. We use the physics convention, i.e., we assume a harmonic time dependence $\u221dexp(\u2212i\omega t)$, where ω denotes the angular frequency. The fields are then expressed as $H(x,t)=\u211c{H\u0302(x)exp(\u2212i\omega t)}$, $E(x,t)=\u211c{E\u0302(x)exp(\u2212i\omega t)}$, where $H\u0302$ and $E\u0302$ denote the complex amplitude vectors. All equations can easily be converted to the engineering convention [expressing the time dependence as $\u221dexp(j\omega t)$] by making the formal substitution $i\u2192\u2212j$. For optical frequencies, typically the permeability is given by its vacuum value μ_{0}. Furthermore, it is here assumed that the dielectric constant can be described by a position dependent complex scalar $\u03f5r(x)$. Under these conditions, Maxwell's equations simplify to

**E**and

**H**. Eqs. (26) and the corresponding boundary conditions define an eigenvalue problem, which yields the electromagnetic resonator modes. General numerical approaches for solving Maxwell's equations include the finite element

^{73}and finite difference time domain

^{35}method, which also have been applied to the simulation of QCL cavities.

^{71,74}

### B. Two-dimensional waveguide model

By eliminating $E\u0302$ from Eqs. (26), a wave equation for $H\u0302$ can be derived.^{75} In many cases, a waveguide geometry is used which does not depend on the longitudinal *x* direction, i.e., $\u03f5r=\u03f5r(y,z)$. Optical outcoupling is then obtained through the cleaved semiconductor facets which serve as partly transparent mirrors. Since the resonator length of up to a few mm is large as compared to the transverse resonator dimensions, the computation of the transverse mode profile in the *yz*-plane can be decoupled from the propagation coordinate and reduces to a 2D problem. The facet transmittance is then calculated based on the obtained transverse field distribution. For a constant waveguide geometry in propagation direction *x*, we can assume a field dependence $Hy,z(x,t)=\u211c{H\u0302y,z(y,z)exp(i\beta x\u2212i\omega t)}$ with the complex propagation constant β, and analogously for the electric field. The wave equation then reduces to two coupled differential equations for the transverse field components

with $k0=\omega \mu 0\u03f50=\omega /c,\u2009p=y,z$, and *q* = *z,y*.^{75} The longitudinal component *H _{x}* is obtained from $\u2207H=0$. Furthermore, the electric field components can be calculated from Eq. (26a). The solution of Eq. (27), along with the corresponding boundary condition of vanishing fields for $y2+z2\u2192\u221e$, constitutes an eigenvalue problem, and the corresponding solutions $H\u0302y,z$ and β correspond to the waveguide modes. As mentioned above, only resonator modes with an electric field component along the

*z*direction are amplified. This is fulfilled for the transverse magnetic (TM) modes, which are characterized by a vanishing magnetic field in the propagation direction, i.e.,

*H*= 0.

_{x}A waveguide mode is frequently characterized using three parameters, the overlap (or field confinement) factor Γ, the waveguide loss coefficient *a*_{w}, and the mirror or outcoupling loss coefficient *a*_{m}, with the total power loss coefficient *a* = *a*_{w} + *a*_{m}. Based on these quantities, the threshold gain is given by $gth=a/\Gamma $. Furthermore, for QCL simulations including the optical cavity field, these parameters enter the simulation of the QCL gain medium to describe the properties of the waveguide mode.^{76,77} Γ, *a*_{w}, and *a*_{m} can be obtained from the mode solutions of Eq. (27). The waveguide loss arises from the absorption in the waveguide layers and is given by $aw=2\u2111{\beta}$. The overlap factor corrects for the fact that the mode only partially overlaps with the gain medium. To reflect the fact that the gain medium only couples to the *E _{z}* component, the overlap factor is defined as

^{78}

where we integrate over the gain medium cross section area *S*_{g} in the enumerator.

Strictly speaking, the calculation of the facet transmission constitutes a full 3D problem, since the facets introduce an abrupt change in *x* direction. However, since the QCL resonator length is large as compared to its transverse dimensions, the computation of the transverse mode profile in the *yz*-plane can be decoupled from the *x* coordinate, as mentioned above. Only for sufficiently wide transverse waveguide dimensions, the facet reflectance *R* can be estimated from Fresnel's formula

with the effective refractive index defined as $neff=\beta /k0$. In general, modal effects lead to an increased reflectance.^{71} Various methods have been developed to extract *R* from the mode solutions provided by Eq. (27).^{79–82} The description of the outcoupling loss by a distributed coefficient *a*_{m} is obtained from $R=exp(\u2212amL)$ with the resonator length *L*, yielding

If one facet is reflection coated and the light is outcoupled only at one side, we obtain $am=\u2212ln(R)/(2L)$.

### C. One-dimensional waveguide slab model

If the waveguide width in lateral *y* direction significantly exceeds its thickness, the waveguide calculations can be reduced to a 1D problem with $\u03f5r=\u03f5r(z)$, corresponding to the simulation of a slab waveguide structure.^{83} This applies for example to typical THz metal-metal waveguide resonators, where the vertical *z* dimension of around 10 *μ*m is often significantly smaller than the lateral *y* dimension $(\u224825\u2212200\u2009\mu m)$.^{71} For TM modes, in addition to *H _{x}* = 0 we have approximately $Hz\u22480$, and the

*y*component is given by $Hy(x,z)=\u211c{H\u0302y(z)exp(i\beta x\u2212i\omega t)}$, where β denotes the complex propagation constant. The one-dimensional wave equation is then obtained from Eqs. (26) or Eq. (27) as

The electric field amplitude is with Eq. (26a) given by

From Eq. (31) we see that both *H _{y}* and $(1/\u03f5r)\u2202zHy$ must be continuous, corresponding to the continuity of the field components

*H*and, with Eq. (32b),

_{y}*E*parallel to the layers.

_{x}Equation (31) can be brought into the same form as the one-dimensional Schrödinger equation (4) by associating $H\u0302y(z),\u2009\u03f5r(z)$, and $\u03f5r\u22121\beta 2\u2212k02$ with $\psi (z),\u2009\u22122m\u2217(z)/\u210f2$, and $E\u2212V(z)$, respectively. Thus, we can use the transfer matrix method introduced in Sec. II D 1 also for solving Eq. (31). We start by dividing the waveguide in *z* direction into layers *n*, $zn\u2264z<zn+\Delta n=zn+1$, with constant relative permittivities $\u03f5r(n)$. Although the gain medium itself consists of different layers, it can be described by its total layer thickness and a single effective dielectric constant since the individual layers are so thin that the electromagnetic wave cannot resolve the structure. For TM modes, the inverse effective dielectric constant is given by $\u03f5r,eff\u22121=(\Delta b\u03f5r,b\u22121+\Delta w\u03f5r,w\u22121)/(\Delta b+\Delta w)$, where $\Delta b$ and $\Delta w$ denote the total thickness of the barriers and wells in the gain medium, respectively, and $\u03f5r,b,\u2009\u03f5r,w$ are the corresponding dielectric constants.^{84,85} In analogy to Eq. (10) for solving the Schrödinger equation, we can write the solution of Eq. (31) in layer *n* as^{83}

with $kn=(\u03f5r(n)k02\u2212\beta 2)1/2$. The propagation through the segment *n* is then described by the matrix $Tn(\Delta n)$ defined in Eq. (14). As mentioned above, Eq. (31) implies the continuity of *H _{y}* and $(1/\u03f5r)\u2202zHy$ across layer boundaries. These matching conditions between two layers can be expressed in terms of matrix Eq. (15) by choosing $k\u0303n=kn/\u03f5r(n)$. There is, however, one fundamental difference between Eqs. (4) and (31): Equation (31) generally has complex eigenvalues $\beta 2$ since $\u03f5r$ becomes complex for materials with optical loss or gain. Consequently, the shooting method described in Sec. II D 1 is not applicable, and a complex root finding algorithm has to be used.

The transmittance *T* = 1 – *R* through the facet can be approximately computed from the mode profile $H\u0302y(z)$. For TM polarization, we obtain from the boundary value method^{79}

with the Fourier transforms

Several extensions of the one-dimensional mode calculations presented above are available. For high aspect ratios of the waveguide cross section, the 2D mode calculation can be approximately reduced to 1D problems in *y* and *z* directions, respectively, by using the effective refractive index method.^{86} A generalization of the transfer matrix scheme for solving the two-dimensional wave equation, Eq. (27), is the film mode matching method, which is especially efficient for waveguides with a rectangular geometry.^{87} Waveguides with a periodic structure in propagation direction, such as a grating for surface emission, can be treated with coupled mode theory.^{34}

## IV. OVERVIEW AND CLASSIFICATION OF CARRIER TRANSPORT MODELS

Depending on the intended goals, carrier transport models with varying degrees of complexity have been used for the QCL gain medium simulation, ranging from simple rate equation approaches to fully quantum mechanical descriptions. The central task is to determine the optical gain, which is proportional to the population inversion between the upper and the lower laser levels, and the current through the heterostructure as a function of the applied bias voltage. Most simulation approaches require the eigenenergies and wave functions of the energy levels in the QCL heterostructure as an input, which are determined by solving the Schrödinger equation or the Schrödinger-Poisson equation system (see Sec. II).

In Fig. 8, various theoretical descriptions of the gain medium with different levels of complexity are illustrated. The most basic model of laser gain media in general is the rate equation approach,^{88} which is also frequently applied to QCLs.^{77,89–94} The electron transport is simply described by transition rates between the relevant energy levels in the gain medium. In the case of the QCL, these are the electron subbands, and the corresponding transitions are referred to as intersubband transitions. The application of the rate equation model to the QCL is illustrated in Fig. 8(a). In addition to the transitions induced by the optical lasing field, nonradiative transitions occur due to various scattering mechanisms, such as the interaction of the electrons with phonons, other electrons or defects in the semiconductor lattice. In the simplest case, these scattering rates are experimental or empirical input parameters.^{89,90} In more advanced implementations, the transition rates are computed based on the corresponding Hamiltonian,^{91,92} only relying on well known material parameters such as the effective mass. Thus, no specific experimental input is required, and no free “fitting” parameters are available. A direct numerical solution is here not possible since the transition rates depend on the *a priori* unknown electron populations. Rather, the steady state solution must be found by simulating the temporal evolution of the system until convergence is reached, or by using an iterative scheme. Such advanced carrier transport methods are generally referred to as self-consistent approaches. Apart from providing an intuitive description, rate equation models are numerically very efficient, and thus are frequently employed for the design and optimization of experimental QCL structures.^{94} An extension are the Maxwell-Bloch equations, which include the carrier-light interaction based on the density matrix formalism.^{95} This approach is used to investigate coherence effects between the laser field and the gain medium, relevant in particular for the pulse formation in mode locked QCLs.^{96–98}

In a quantum well heterostructure such as the QCL gain medium, the electrons are only confined in the growth direction. Thus, energy quantization occurs in one dimension, while the electrons can still freely move in two dimensions. Due to this free in-plane motion in the quantum well, the electrons have kinetic energy in addition to the eigenenergy of the corresponding quantized level, as illustrated in Fig. 8(b). For a more detailed modeling and an improved understanding of the carrier transport in QCLs, the in-plane motion of the electrons must be taken into account. Such three-dimensional (3D) simulation approaches also consider intrasubband transitions occurring between different kinetic energies within a subband, and yield the electron distribution in each subband in addition to the level populations. On the other hand, the complexity increases significantly as compared to one-dimensional (1D) descriptions such as above discussed rate equation approach, since transitions are now characterized by the initial and final subbands as well as the corresponding kinetic energies. An example for a self-consistent 3D approach is the EMC method. Here, the intrasubband processes are fully taken into account, and the scattering rates are self-consistently evaluated based on the corresponding Hamiltonian. Combining versatility and reliability with relative computational efficiency, the EMC approach has been widely used for the analysis, design, and optimization of QCLs.^{46,47,99–110}

In both the rate equation and EMC approach, the carrier transport is described by scattering-induced transitions of discrete electrons between the quantized energy levels, corresponding to a hopping transport model. Such methods are referred to as semiclassical, since the energy quantization in the QCL heterostructure is considered, but quantum coherence effects and quantum mechanical dephasing are not included. Various quantum transport simulation approaches have been developed which take into account quantum correlations between the energy levels. One example is the density matrix method. Its 1D version can be seen as a generalization of the rate equation approach, and is frequently used for the analysis and optimization of THz QCLs.^{17,111–113} In its 3D form^{114,115} illustrated in Fig. 8(c), it corresponds to a generalization of semiclassical approaches based on the Boltzmann transport equation such as EMC.^{114} Also the nonequilibrium Green's function (NEGF) method, considered the most general quantum transport approach, has been applied to the simulation of QCL structures.^{116–124} Quantum transport approaches are numerically much more demanding than their semiclassical counterparts. On the other hand, quantum coherence effects can play a pronounced role especially in THz QCLs where the energetic spacing between the quantized levels is relative small,^{115,125} while they are less relevant in mid-infrared QCLs.^{114}

In Table I, a classification of the different simulation approaches covered in this review is given. Here, we differentiate between semiclassical schemes based on hopping transport between the quantized energy levels, and quantum transport approaches taking into account quantum correlations. We divide the methods into empirical approaches relying on empirical or experimental input parameters, and self-consistent schemes which evaluate the transition rates based on the corresponding Hamiltonian. Furthermore, we distinguish between 1D modeling techniques only considering the subband populations and intersubband transitions, and 3D approaches also taking into account the electron distribution in the subbands and the intrasubband dynamics.

. | Semiclassical . | Quantum transport . |
---|---|---|

Empirical | ||

1D | Rate equations (V A) | 1D density matrix (VII) |

Maxwell-Bloch ^{a} (V B) | ||

Self-consistent | ||

1D | Rate equations (VI E) | 1D density matrix (VII) |

3D | Monte Carlo (VI F) | 3D density matrix (VII) |

NEGF (VIII) |

. | Semiclassical . | Quantum transport . |
---|---|---|

Empirical | ||

1D | Rate equations (V A) | 1D density matrix (VII) |

Maxwell-Bloch ^{a} (V B) | ||

Self-consistent | ||

1D | Rate equations (VI E) | 1D density matrix (VII) |

3D | Monte Carlo (VI F) | 3D density matrix (VII) |

NEGF (VIII) |

^{a}

Only the carrier-light interaction is modeled using a density matrix formalism, while scattering is treated based on rate equations.

## V. EMPIRICAL APPROACHES

The most basic approach for modeling the electron dynamics in a laser is to use experimental or empirical transition rates between the relevant energy levels in the laser gain medium.^{88} For QCLs, these levels correspond to the quantized eigenstates of the heterostructure, which have to be found by solving the Schrödinger equation, Eq. (4). The electron dynamics is then described by rate equations.^{88} Often, only the nonradiative transitions are considered which occur due to various scattering mechanisms, such as the interaction of the electrons with phonons, other electrons or defects in the semiconductor lattice. Simulations not including the optical cavity field can be used to investigate under which conditions sufficient optical gain is obtained, so that lasing operation can start at all. In this way, parameters such as the threshold current density and maximum operating temperature can be extracted. To investigate the lasing operation itself, including optical output powers, saturation effects or the intrinsic linewidth,^{90} the optical field has to be included as well. Modified empirical scattering-rate approaches have been developed for specific types of QCLs.^{126} Maxwell-Bloch equations which are a generalization of the rate equation approach can be used to model the coherent interaction between the laser field and the gain medium,^{95} e.g., to investigate the formation of optical instabilities in mode-locked QCLs.^{96}

### A. Empirical rate equations

The rate equations for a laser are given by^{88}

The first two terms contain the relaxation transitions due to scattering, e.g., the interaction of the electrons with phonons, other electrons or defects in the semiconductor lattice. Also spontaneous emission can be included here. The scattering rate from a level *j* to *i* is often expressed in terms of an inverse scattering lifetime $\tau ji\u22121$, where $\tau i\u22121=\u2211j\u2260i\tau ij\u22121$ indicates the total inverse lifetime of level *i*. The last sum contains the lasing transitions, where $Wijopt$ are the stimulated optical transition rates for those transitions where an optical field at or near the corresponding frequency $|\omega ij|=|Ei\u2212Ej|/\u210f$ is present.^{88} The rates $Wijopt$ are proportional to the optical intensities in the corresponding lasing modes. Typically, only one or a few transitions contribute to lasing. Furthermore, $nis$ is the electron sheet density of subband *i*, i.e., the electron number divided by the in-plane cross section area *S*. This quantity is often used in QCL heterostructures where energy quantization occurs in one dimension and the electrons can still freely move in in-plane direction.

Commonly, the QCL heterostructure is designed strictly periodically, as illustrated in Fig. 1. Apart from fabrication tolerances, a periodic model is valid for the central QCL periods far away from the contacts, if effects such as domain formation^{127} and local variations of the optical field intensity can be neglected. The sum in Eq. (36) can then be restricted to one representative central period

with $i,j=1..N$, where *N* is the number of subbands in each period. Here, $\tau \u0302ji\u22121=\u2211n\u2208\mathbb{Z}\tau j,i+nN\u22121$ includes the transitions to all equivalent levels in the different periods, and analogously for $W\u0302jiopt$. The total sheet density in each period is determined by the doping sheet density

The steady state solution is obtained by setting $dtnis=0$. If we are not interested in the lasing operation itself, but only if sufficient inversion for lasing is obtained, stimulated optical effects can be excluded, $W\u0302ijopt=W\u0302jiopt=0$. The subband populations $nis\u2009(i=1..N)$ can then be found by solving the linear equation system Eq. (37) with $dtnis=0$ and $i=1..(N\u22121)$, complemented by Eq. (38) to obtain a linearly independent system and thus a unique solution. To include lasing, this system has to be complemented by equations describing the optical intensities in the lasing modes, since the $W\u0302ijopt$ is intensity dependent.

Often the empirical rate equation is restricted to the most important subbands and transitions to obtain compact analytical results. A common model is the three-level system,^{89} illustrated in Fig. 9. This description only includes the upper laser level 2, lower laser level 1, and a reservoir level 0 representing the extraction and injector subbands. A fraction $\eta 2$ of the current density *J* is injected into the upper laser level, and $\eta 1$ into the lower laser level. Furthermore, only the transition rates $\tau 21\u22121,\u2009\tau 20\u22121$, and $\tau 10\u22121$ are included, with the inverse upper and lower laser level lifetimes $\tau 2\u22121=\tau 21\u22121+\tau 20\u22121,\u2009\tau 1\u22121=\tau 10\u22121$. With Eq. (36), the rate equations for the upper and lower laser level are then obtained as

where $\sigma 21opt$ is the cross section for stimulated emission from level 2 to 1 and *I* is the optical intensity of the laser radiation at frequency $\omega 21=(E2\u2212E1)/\u210f$. If lasing is neglected, *I* = 0, the steady state population inversion is with d_{t} = 0 in Eq. (39) obtained as^{89}

The three-level model can be extended to take into account further effects such as thermal backfilling of the lower laser level and backscattering from level 1 to level 2 in terahertz QCLs, which can be modeled by introducing additional lifetimes $\tau 01$ and $\tau 12$, respectively.^{90}

### B. Maxwell-Bloch equations

Bloch equations are a generalization to the rate equation approach where the interaction of the laser level electrons with the optical field is modeled using a density matrix formalism rather than scattering rates.^{95} In this way, optical nonlinearities and coherence effects between the laser field and the gain medium can be considered, while the carrier transport due to nonradiative mechanisms is included by the corresponding lifetimes as for the rate equations. To describe the optical propagation, the Bloch equations are complemented by Maxwell's or related equations, such as the wave equation. This model has, for example, been used to investigate the formation of optical instabilities in QCLs,^{96,97} and to study the possibility of self-induced transparency modelocking.^{98} The nonlinear optical dynamics in the gain medium also plays an important role for the recently demonstrated QCL-based frequency combs.^{128}

The interaction of a classical optical field with a two-level system is in the density matrix formalism described by^{95}

where the asterisk denotes the complex conjugate. $Ez(x,t)$ represents the field component in the growth direction *z* of the heterostructure since the other components do not interact with the gain medium, and $d12=\u2212e\u27e81|z|2\u27e9$ is the corresponding dipole matrix element of the laser transition. Furthermore, $\rho ij(x,t)=\rho ji\u2217(x,t)=\u27e8i|\rho \u0302(x,t)|j\u27e9$ are the density matrix elements, and $\Delta 21(x,t)=\rho 22(x,t)\u2212\rho 11(x,t)$ is the inversion. The density matrix can be normalized so that $\rho ii=nis/ns$ gives the relative population of subband *i*. Dissipative processes are phenomenologically included by adding decay terms with relaxation rates $\gamma E$ and $\gamma 21$, describing the energy relaxation and dephasing, respectively. $\Delta 21eq$ is the equilibrium inversion which the system approaches for *E _{z}* = 0. If the gain medium is not homogeneous along the propagation direction, as is the case for self-induced transparency mode locking structures,

^{98}the parameters

*d*

_{12}, $\gamma E,\u2009\gamma 21,\u2009\omega 21$, and $\Delta 21eq$ depend on the propagation coordinate

*x*.

Assuming weak nonlinearity and inhomogeneity, the optical field propagation can be described by the wave equation.^{95} For propagation in *x* direction, it is given by

$Pz(x,t)$ is the polarization component in *z* direction due to the lasing transition, given by $Pz=(ns/Lp)(d12\rho 21+d12\u2217\rho 21\u2217)$. Here, *L*_{p} is the length of a single QCL period and *n*^{s} /*L*_{p} thus corresponds to the average electron concentration.^{97} Furthermore, *n*_{0} denotes the refractive index of the gain medium material.

*E _{z}* and $\rho 21$ are typically expressed by their slowly varying envelope functions

where $k21=\omega 21n0/c$ and c.c. denotes the complex conjugate. Inserting Eq. (43) into Eq. (41) and neglecting the rapidly oscillating terms $\u221dexp(\xb12i\omega 21t)$, we obtain the density matrix equations in the rotating wave approximation^{95}^{,}

^{95,97}

furthermore assuming $|\u2202t2\eta 21|,|\omega 21\u2202t\eta 21|\u226a|\omega 212\eta 21|$. Here, the overlap factor Γ, Eq. (28), has been added to correct for the fact that the optical mode only partially overlaps with the gain medium, and a loss term with the power loss coefficient *a* has been included. For applying the Maxwell-Bloch equation model to QCLs, typically the simplifying assumption is made that the lower laser level is depopulated very efficiently, $\tau 1\u21920$, resulting in $n1s=0$ in Eq. (39).^{97} The inversion is then directly given by the upper laser level population, $\Delta 21=n2s$. By comparison of Eq. (44) with Eq. (39), we obtain $\Delta 21eq=\tau 2\eta 2J/(nse),\u2009\gamma E=\tau 2\u22121$.

The Maxwell-Bloch equations, Eqs. (44) and (45), are a versatile approach to describe the optical dynamics in QCLs, and additional effects can straightforwardly be implemented. For example, the influence of spatial hole burning due to the standing wave modes in a linear cavity has been extensively studied,^{96,97,129} and dispersion as well as saturable absorption have been added to the model.^{130}

#### 1. Optical gain and transition rates

The Bloch equations Eq. (44) can be used to derive the optical gain coefficient and transition rate associated with stimulated emission and absorption.^{95} For a monochromatic electromagnetic field, the stationary solution of Eq. (44a) is obtained by setting $\u2202t\eta 21=0$

Multiplying Eq. (45) from left with $E\u0302z\u2217$ and adding the complex conjugate, we obtain with Eq. (46)

with the power gain coefficient

Here, we have replaced the electric field by the optical intensity $I=\u03f50cn0|E\u0302z|2/2$. The optical power inside the resonator is given by $P=ISg/\Gamma $, where *S*_{g} is the cross section area of the gain medium and $Sg/\Gamma $ corresponds to the effective area of the waveguide mode. Frequently, *I* is assumed to change only slightly along the resonator, i.e., the intensity is averaged over the *x* coordinate. This assumption is valid especially for the case of moderate output coupling at the facets where the mirror loss can be described by a distributed coefficient *a*_{m}, Eq. (30). Equation (47) then simplifies to

The transition rate due to the optical field is with $\u2202t\Delta 21|opt=2\u2202t\rho 22|opt=\u22122\u2202t\rho 11|opt$ obtained from Eqs. (44b) and (46) as

The contribution $\u221d\rho 22I$ leading to a reduction of $\rho 22$ is due to stimulated emission, while the contribution $\u221d\rho 11I$ corresponds to absorption. For a slightly detuned optical field at a frequency $\omega \u2260\omega 21$, Eqs. (48) and (50) can be adapted by replacing $\gamma 21\u22121$ with $\pi L(\omega )$, where $L(\omega )$ is the Lorentzian lineshape function^{95}

Thus, we obtain with the sheet densities $n1s=ns\rho 11$, $n2s=ns\rho 22$,

and

Comparison with Eq. (36) yields for the stimulated optical transition rates

## VI. SELF-CONSISTENT SEMICLASSICAL APPROACHES

Advanced self-consistent simulation approaches only rely on well known material parameters such as the effective mass, and no further specific experimental or empirical input is required. This also means that no adjustable parameters are available to fit the simulation results to experimental data. These approaches are based on the evaluation of the transitions between the various states due to different scattering mechanisms, including the interaction of the electrons with phonons, impurities, and other electrons. The associated scattering rates are computed based on the corresponding Hamiltonian. A direct numerical solution of the resulting equations is not possible since the transition rates depend on the initially unknown electron populations. Rather, the steady state solution must be self-consistently found by simulating the temporal evolution of the system until convergence is reached, or by using iterative schemes. These methods rely on the subband wave functions and eigenenergies found by solving the Schrödinger or Schrödinger-Poisson equation, as described in Sec. II. The carrier transport is then modeled by transitions of discrete electrons between the quantized energy levels, also referred to as hopping transport. Thus, these methods are called semiclassical, since the energy quantization in the QCL heterostructure is considered, but quantum coherence effects and dephasing mechanisms are not included. Formally, semiclassical carrier transport descriptions can be derived from the more general density matrix formalism by neglecting the contribution of the off-diagonal matrix elements, i.e., only considering the diagonal elements corresponding to the occupations of the states.^{114,131} In the following, the most relevant scattering mechanisms in the QCL will be discussed. Furthermore, the self-consistent rate equation approach and the EMC method will be described, which are the two most widely used advanced semiclassical QCL simulation schemes.

### A. Scattering mechanisms and transition rates

A transition of a carrier from one state to another due to a perturbation is referred to as a scattering process. The perturbation is described by a corresponding potential *V*, which can be static or time dependent.^{61} This results in different classes of scattering processes, illustrated in Fig. 10 for the case of one-dimensional electron confinement as in QCL heterostructures. Elastic scattering, shown in Fig. 10(a), occurs for time constant potentials. Here, the carrier energy is conserved. Relevant mechanisms in QCLs include impurity, interface roughness, and alloy scattering. For potentials with harmonic time dependence, $V\u221dcos(\omega 0t)$, the carrier energy is changed by $\u2213\u210f\omega 0$, corresponding to the case of emission and absorption, respectively. This is referred to as inelastic scattering [see Fig. 10(b)]. An important inelastic mechanism is optical phonon scattering, and also the interaction with photons can be viewed as an inelastic scattering process. A special case is intercarrier process illustrated in Fig. 10(c) such as electron-electron scattering, where two electrons are involved in the scattering event. Figure 11 shows the influence of various scattering mechanisms on the spectral gain obtained from an EMC simulation for two different terahertz QCL designs.^{107}

QCLs are n-type devices which are typically based on direct bandgap semiconductors. Thus, in the following, we restrict our treatment of scattering to the conduction band Γ valley where also the lasing transitions occur. However, we note that under certain conditions transitions to other valleys can affect QCL operation.^{45–48} Scattering causes an electron transition from an initial state $|ik\u27e9$ to a final state $|jk\u2032\u27e9$ in the QCL heterostructure, where $k=(kx,ky)T$ and $k\u2032=(kx\u2032,ky\u2032)T$ are the corresponding in-plane wave vectors. The states are described by their wave functions Eq. (2) and eigenenergies Eq. (5). For the initial state, we have $\psi 3D,i=S\u22121/2\psi i(z)exp(ikr)$ with $r=(x,y)T$ and $Eik=Ei+\u210f2k2/(2mi\u2225)$ where $\psi i(z)$ and *E _{i}* are obtained from Eq. (4) as described in Sec. II A, and analogously for the final state.

For elastic scattering processes, *V* is constant. The corresponding matrix element is defined as

For inelastic processes, we assume a harmonic potential of the form $V=V0exp(iQx\u2212i\omega 0t)+V0\u2217exp(\u2212iQx+i\omega 0t)$, where **Q** is, for example, the phonon wave vector. The matrix element can in analogy with Eq. (55) be obtained as

Assuming bound wave functions $\psi i,j(z)$ in Eqs. (55) and (56), the integration over the *z* coordinate can be taken from $\u2212\u221e$ to $\u221e$. This is also consistent with treating the gain medium as an infinitely extended periodic heterostructure. Furthermore, the in-plane cross section *S* is assumed to be macroscopic, and thus the integration can for scattering rate calculations be extended from $\u2212\u221e$ to $\u221e$ also in *x* and *y* direction.

The transition rate from an initial state $|ik\u27e9$ to a final state $|jk\u2032\u27e9$ is obtained from Fermi's golden rule,^{61} given by

for elastic scattering processes and

for inelastic processes. Here, $Wik,jk\u2032+$ corresponds to the emission rate, caused by the component $exp(i\omega 0t)$, and $Wik,jk\u2032\u2212$ refers to absorption due to the component $exp(\u2212i\omega 0t)$. The Dirac δ function ensures energy conservation. For elastic processes, we obtain from $Ejk\u2032=Eik$

Analogously, energy conservation yields for inelastic scattering

The computation of the total transition rates from an initial state $|ik\u27e9$ to a final subband *j* involves the summation over wave vectors. These sums can be converted to integrals, introducing a factor of $Ld/(2\pi )$ per dimension where the device length *L _{d}* in the corresponding direction is assumed to be large enough that quantization effects can be neglected.

^{132}An example is the summation over the final in-plane wave vector $k\u2032$ which is two-dimensional. It can furthermore be advantageous to express $k\u2032$ in polar coordinates $|k\u2032|$ and $\varphi $, and introduce a kinetic energy variable $Ejkin=\u210f2|k\u2032|2/(2mj\u2225)$. Thus, we obtain

Spin degeneracy is not considered here, since for single-electron scattering processes the spin is conserved. For electron-electron scattering, the spin degeneracy must, however, be taken into account, as more closely discussed in the corresponding section. Analogously, summation over a three-dimensional wave vector, such as the phonon wave vector **Q**, can in a crystal lattice of volume $\Omega c$ be approximated by

### B. Phonon scattering

A phonon is a quasiparticle associated with the lattice vibrations in a crystal, representing an excited quantum mechanical state in the quantization of the vibrational modes. Classically speaking, for an atom located at position **x**, lattice vibrations are described by a displacement vector $u=Usin(Qx\u2212\omega Qt)$, with the amplitude **U**, wave vector **Q**, and angular frequency $\omega Q$. The normal modes are the solutions of **u** for which the lattice uniformly oscillates at a single frequency $\omega Q$, and a phonon corresponds to an elementary vibrational motion. The associated relation between wave vector **Q** and frequency $\omega Q$, the so-called dispersion relation, defines a phonon branch, as illustrated in Fig. 12. Acoustic modes are sound waves where two consecutive atoms move in the same direction, and we have $\omega Q=0$ for **Q = 0**. For optical modes, two consecutive atoms in the same unit cell move in opposite direction, and $\omega Q$ at **Q = 0** typically corresponds to infrared optical frequencies. Furthermore, the wave propagation can be predominantly longitudinal $(Q\u2225U)$ or transverse $(Q\u22a5U)$; the corresponding branches are then referred to as longitudinal or transverse branches. In three dimensions, there is a single longitudinal acoustic (LA) and 2 transverse acoustic (TA) branches as well as *N*_{u} – 1 longitudinal optical (LO) and $2(Nu\u22121)$ transverse optical (TO) branches, where *N*_{u} denotes the number of atoms per unit cell. For example, GaAs has two atoms per unit cell, i.e., *N*_{u} = 2.

The lattice vibrations lead to a perturbation of the carriers and thus carrier scattering. Depending on the mechanisms, there are different types of phonon scattering. Non-polar phonon scattering occurs in all crystals and is due to acoustic and (for $Nu\u22652$) TO phonons. Here, the lattice vibrations lead to a time dependent change of the conduction (and valence) band energy. Polar scattering only occurs in polar semiconductors, e.g., III–V semiconductors such as GaAs. It is due to LO phonons, where the out-of-phase movement of the neighbouring atoms of different types causes a local dipole moment, resulting in oscillating electric fields. Covalent semiconductors such as group IV materials do not exhibit polar scattering.

The dominant phonon scattering mechanism in QCLs is due to LO phonons. Since the optical phonon deformation potential approaches zero at the conduction band Γ point due to spherical symmetry, TO phonon scattering is negligible for the Γ valley; it can however play a significant role for intervalley transitions and intravalley scattering in the *X* and valence band Γ minimum.^{133–135} Also acoustic phonons tend to play a secondary role in QCLs.^{136}

#### 1. Non-polar phonon scattering

Non-polar phonon scattering occurs in all crystals and is due to acoustic and (for $Nu\u22652$) TO phonons. The energy of the vibrating atom is the sum of its potential and kinetic energy, which is equal to twice its average kinetic energy. Thus, the total energy of the vibrational mode becomes for $u=Usin(Qx\u2212\omega Qt)$

where the number of atoms with mass *m* in a crystal lattice of volume $\Omega c$ and density $\rho c$ is $Na=\Omega c\rho c/m$. Thus, we obtain the amplitude for a mode occupied by a single phonon of energy $E=\u210f\omega Q$

##### A. Acoustic phonons

Since phonon confinement does not play an important role for the cascade structures (see discussion in Sec. VI B 2), scattering from bulk phonons obeying Bose distribution is considered. Lawaetz has shown in Ref. 137 that the deformation potential method of Bardeen and Shockley^{138} may be applied to the scattering of electrons with acoustic phonons. This method is basically a Taylor expansion of the scattering potential in the phonon momentum **Q**. In the case of vanishing screening, Kittel and Fong derive in Ref. 139 the change of the electronic energy in lowest order of the lattice deformation. Since for acoustic phonons, adjacent atoms move in the same direction, regions of compression or dilatation extend over several lattice sites, and the crystal can be described as an elastic continuum. The resulting strain gives rise to a time dependent change of the conduction band energy $V=\delta Vc$. The corresponding constant of proportionality is the deformation potential Ξ (for valleys without rotational symmetry, Ξ becomes a tensor), i.e., $V=\Xi \u2207u$ for small displacements.^{61} We obtain

i.e., only LA phonons contribute since **QU** = 0 for transverse phonons. For acoustic phonons, we have $\omega Q=0$ for $Q=|Q|=0$, and the dispersion relation can be described by the linear approximation $\omega Q=vsQ$ for small *Q*, where *v*_{s} is the longitudinal sound velocity (see Fig. 12). With Eq. (64), the perturbation potential is thus

with the amplitude

Equation (58) gives the transition rate for an electron in an initial state $|ik\u27e9$, which is scattered to a final state *j* with a wave vector $k\u2032$ by a phonon with wave vector $Q=[q,qz]T$, where **q** corresponds to the in-plane component. Using the definition of the three-dimensional wave functions in Eq. (2), the corresponding matrix element of the perturbation potential Eq. (66) is with Eq. (56) given by

with the form factor

For bound states, $\psi i,j$ can be chosen to be real, and we have $|Fji+(qz)|2=|Fji\u2212(qz)|2=|Fji(qz)|2$. The scattering rate obtained from Eq. (58) is thus

where $Wik,jk\u2032+$ and $Wik,jk\u2032\u2212$ refer to emission and absorption, respectively. Here, we have used that $|\u222bSexp[i(k\xb1q\u2212k\u2032)r]d2r|2$ can be approximated by $4\pi 2S\delta (k\u2213q\u2212k\u2032)$ for sufficiently large in-plane cross sections *S*. Since *M*(*Q*) refers to a single phonon, the phonon occupation number of a mode in thermal equilibrium has been added in Eq. (70), given by the Bose-Einstein distribution

For emission, the factor *N*_{Q} + 1 is used to also include spontaneous emission processes.

The total transition rate from a given initial state $|ik\u27e9$ to a subband *j* is obtained by summing over all wave vectors $k\u2032$ and **Q**. These sums can be converted to integrals using Eqs. (61) and (62). With Eqs. (70) and (67), we thus obtain the total transition rate from a given initial state $|ik\u27e9$ to a subband *j*^{132}

Since acoustic phonons do not carry much energy, often the quasi-elastic approximation is applied, treating acoustic phonon scattering as elastic process. This is achieved by neglecting the phonon energy term $\xb1\u210f\omega Q$ in Eqs. (70) and (72). Besides, we can approximate Eq. (71) as

for all but the lowest temperatures, which is referred to as equipartition approximation. Thus, we obtain

With Eqs. (61) and (62), we furthermore obtain $Wik,j\xb1(Q)=0$ for $Ej>Eik$, and otherwise

Thus, the total transition rate is given by

which includes both emission and absorption and thus contains an additional factor of 2.

##### B. Transverse optical phonons

Transverse optical phonons are treated in a similar manner as acoustic phonons. Here, adjacent atoms move in opposite directions. Thus, the crystal vibrations cannot be described as an elastic continuum as for acoustic phonons, but correspond rather to oscillations of the two sublattices with respect to each other. As a consequence, the potential change is directly a function of the displacement vector **u**, $\delta Vc=DTOu=DTOu$. Here, **D**_{TO} is the optical deformation potential field which assumes a different value for different kinds of intra- and intervalley transitions, and *D*_{TO} is the component parallel to **u**. For optical phonons, $\omega Q$ reaches an extremum $\omega Q=\omega TO$ for *Q* = 0, and for small *Q*, the dispersion relation can thus be approximated by $\omega Q=\omega TO$ (see Fig. 12). The scattering rate for an electron from an initial state $|ik\u27e9$ to a final state $|jk\u2032\u27e9$ due to phonons with wave vector **Q** is again given by Eq. (70), with

In analogy to Eq. (72), we obtain the total transition rate

Approximating $\omega Q=\omega TO$, the phonon number *N*** _{Q}** in Eq. (79) does not depend on

**Q**, and thus we can set

*N*

**=**

_{Q}*N*. With Eq. (76), the result then simplifies to

_{Ph}for $Ej\u2264Eik\u2213\u210f\omega TO$, and becomes 0 otherwise.

#### 2. Longitudinal optical phonons

As for TO phonons, adjacent atoms move in opposite directions for LO phonons. Again the dispersion relation can be approximated with a constant value $\omega Q=\omega LO$ for small *Q* (see Fig. 12), where $\u210f\omega LO$ is the LO phonon energy at *Q* = 0. Polar semiconductors contain different types of atoms carrying effective charges. For two different atoms with masses *m*_{A} and *m*_{B} in a unit cell, these effective charges are given by $\xb1qeff$ with $qeff2=\u03f50\omega LO2\mu nc\u22121(\u03f5r,\u221e\u22121\u2212\u03f5r,0\u22121)$.^{61} Here, $\mu \u22121=mA\u22121+mB\u22121$ is the inverse reduced mass, and *n*_{c} is the number of cells per unit volume. Furthermore, $\u03f5r,0$ denotes the static dielectric constant and $\u03f5r,\u221e$ is the dielectric constant at very high frequencies where the ions cannot respond to a field anymore and thus only the electronic polarization remains. For LO waves, the out-of-phase movement of the atoms causes a local dipole moment $p(t)=qeffu(t)$, where $u=Usin(Qx\u2212\omega Qt)$ now indicates the relative displacement of the atoms in the unit cell; by contrast, transverse waves do not generate dipoles. The amplitude for excitation by a single phonon is given by Eq. (64) where $\rho c$ is replaced by μ*n*_{c}. The resulting ionic polarization is $Pion=pnc$, and the electric field **E** can be obtained from the material equation $D=\u03f50E+Pion$ with **D = 0** due to the absence of free carriers. With the potential $V=e\u222bEdx$, the resulting scattering rate from a state $|ik\u27e9$ to $|jk\u2032\u27e9$ is again given by Eq. (70), where

The total transition rate from a given initial state $|ik\u27e9$ to a subband *j* is again obtained by summing over all wave vectors $k\u2032$ and **Q**, using Eqs. (61) and (62). For $Ej\u2264Eik\u2213\u210f\omega LO$, we obtain

and $Wik,j\xb1=0$ otherwise. Here,

with the in-plane phonon wave vector magnitude

where $k\u2032$ is given by Eq. (60) with $\omega 0=\omega LO$.

##### A. Corrections due to screening, hot phonons, phonon modes, and lattice heating

In a heterostructure, the assumption of bulk phonon modes is justified for acoustic phonons, but not necessarily for LO phonons. The LO phonon modes arising from different dielectric constants in the well and barrier material can be described by microscopic and macroscopic approaches. The macroscopic dielectric continuum model, which is widely used in this context, yields slab modes which are confined in each layer, and interface modes corresponding to exponential solutions around the interfaces.^{140} However, it has been shown that the bulk mode approximation is valid for most QCL designs.^{48,101} Specifically, for not too different compositions of the barrier and well material and thus similar values of $\u03f5r,0$, the use of bulk modes is a good approximation.^{101}

Screening can be included by changing Eq. (81) to^{141,142}

where *q*_{s} is the inverse screening length. For Eq. (83), we then obtain with $q\u0303=q2+qs2$

The simplest model is to use Debye screening with the bulk inverse Debye screening length

where *n*_{e} is the averaged three-dimensional electron density. However, Debye screening is generally only valid for high temperatures,^{143} since *q*_{s} diverges for $TL\u21920$.

Deviations from the equilibrium phonon distribution Eq. (71) due to phonon emission and absorption in the heterostructure, also referred to as hot phonons, are frequently considered in simulations.^{104,144} The hot phonon effect is commonly incorporated by using a relaxation time approach, describing the decay of LO phonons into other phonon modes by a corresponding lifetime.^{145}

Typically, in simulations it is assumed that the lattice temperature of the gain medium is identical to a given heat sink temperature or to the room temperature. However, this approximation is only legitimate for pulsed operation at low duty cycle, and generally fails for continuous wave operation.^{146} Lattice heating due to the dissipated electrical power can be self-consistently modeled by coupled carrier transport and thermal simulations,^{147,148} yielding the actual temperature profile. The temperature distribution is computed by solving the heat diffusion equation, where the heat generation rate is obtained from the carrier transport simulation based on the dissipated electrical power or the LO phonon scattering rate.

### C. Electron-electron scattering

Electron-electron scattering arises from the collision of electrons with other electrons. More precisely, this scattering mechanism can be divided into two main contributions, the short-range interaction between two electrons via the screened Coulomb potential and the long-range electron-plasmon coupling.^{149} The latter is usually neglected in QCL simulations and is also not considered here, since it only becomes relevant at elevated doping levels.

Many-electron effects are to lowest order already considered in the Poisson equation, Eq. (20), corresponding to a mean field treatment of charges. In the semiclassical treatment, an inclusion of higher order effects is obtained by corresponding scattering rates. Electron-electron scattering is much more computationally demanding than the single-electron processes, hampering its inclusion in quantum mechanical simulations of QCLs beyond the mean-field approximation.^{118}

In semiclassical simulations, electron-electron scattering is typically implemented as a two-electron process.^{141,150} An electron in an initial state $|ik\u27e9$ scatters to a final state $|jk\u2032\u27e9$, accompanied by a transition of a second electron from a state $|i0k0\u27e9$ to $|j0k0\u2032\u27e9$. Nonparabolicity effects have for example been considered by Bonno and Thobel.^{55} To facilitate the treatment, we assume here an identical mass *m** for all subbands. The total scattering rate from $|ik\u27e9$ to a subband *j* is then obtained as^{141,150}

with the in-plane cross section area *S* and carrier distribution function $fi0(k0)$ for the state $|i0k0\u27e9$. θ is the angle between $g=k0\u2212k$ and $g\u2032=k0\u2032\u2212k\u2032$, and $q=k\u2212k\u2032$ (with $q=|q|$) denotes the exchanged wave vector. For the definition of the Coulomb matrix element, different conventions exist.^{141,151} Here, we use the widely employed convention of Goodnick and Lugli,^{110,141,150,152} where the bare Coulomb matrix element is with $\u03f5=\u03f50\u03f5r,0$ given by

Electron-electron scattering is sensitive with respect to the spin. Furthermore, screening has to be considered, caused by the response of the other electrons to changes in the electrostatic field associated with the scattering process. Often, only scattering of electron pairs with antiparallel spin is considered, since the contribution of carriers with parallel spin is smaller.^{141,153} Furthermore, basic models consider screening by a constant screening wave number *q*_{s}. For the transition matrix element $Mii0jj0$, we then obtain

where the factor $q/(q+qs)$ corrects the bare Coulomb matrix element for screening effects, and the factor 1/2 arises from neglecting parallel spin contributions. From energy conservation, we obtain for the exchanged wave vector with $g=|g|$ and $g02=4m\u2217(Ei+Ei0\u2212Ej\u2212Ej0)/\u210f2$

A frequently used basic screening model is the single subband screening approximation, assuming that most electrons in each period reside in a single subband *i*, corresponding to the ground state. The obtained screening wave number is then^{110,152}

#### 1. Advanced spin and screening models

Advanced implementations of spin and screening effects have been developed for a more accurate description of electron-electron scattering.^{68} Screening models with varying degrees of sophistication are employed to obtain the screened Coulomb matrix elements $Vii0jj0s$ from $Vii0jj0b$ in Eq. (89). In the random phase approximation (RPA), $Vii0jj0s(q)$ is found by solving the equation system^{154}

Here, $\Pi mn(q)$ is the polarizability tensor, given in the long wavelength limit $(q\u21920)$ by

Commonly, simplified screening models using a constant screening wave number *q*_{s} are employed to avoid the numerical load associated with solving Eq. (93).^{110,152}$Vii0jj0s$ is then obtained by replacing the prefactor $e2/(2\u03f5q)$ in Eq. (89) by $e2/[2\u03f5(q+qs)]$. This implementation is also used in Eq. (90). The single subband screening model, Eq. (92), can be obtained from Eq. (93) by assuming that most of the electrons reside in a single subband, and screening is caused only by this subband.^{110,152} A somewhat improved approach which considers all subbands equally is the modified single subband model^{152} with

where *i* sums over the subbands in one period. An alternative approach consists in treating intersubband scattering as unscreened, thus applying Eq. (95) only to the intrasubband matrix elements.^{46}

For collisions of electrons with parallel spin, interference occurs between $Vii0jj0s$and the “exchange” matrix element $Vii0j0js$.^{153} Accounting for this exchange effect, we obtain^{150,153}

where

and *p*_{a} = *p*_{p} = 1/2 are the probabilities for antiparallel and parallel spin collisions, respectively. There are two common approaches to implement electron-electron scattering without explicitly considering the spin dependence. One method which has been used in Eq. (90) is to completely neglect the parallel spin collisions,^{141,153} implying *p*_{a} = 1/2, *p*_{p} = 0 in Eq. (96). This approach tends to overestimate the exchange effect. Another common method is to ignore the exchange effect, i.e., parallel spin collisions are treated the same way as antiparallel spin contributions.^{153} This corresponds to *p*_{a} = 1, *p*_{p} = 0 in Eq. (96).

Figure 13 contains EMC simulations of the gain spectra for different implementations of screening and spin effects.^{68} In Fig. 13(a), the exchange effect is fully included, but different screening models are employed. The reference curve is based on the exact evaluation of the RPA (solid curve). Applying the simplified screening model Eq. (95) to all matrix elements overestimates the screening of the intersubband elements, thus leading to an underestimation of scattering. The resulting spectral gain profile (dashed curve) features an excessively narrow gain bandwidth and enhanced peak gain. On the other hand, completely ignoring the screening effect for the intersubband matrix elements overestimates the intersubband scattering, thus resulting in a lowered and broadened gain profile (dotted curve). In Fig. 13(b), screening is included in the RPA. Results are shown for *p*_{a} = *p*_{p} = 1/2 (solid curve), *p*_{a} = 1, *p*_{p} = 0 (dotted curve), *p*_{a} = 1/2, *p*_{p} = 0 (dashed curve), and *p*_{a} = *p*_{p} = 0 (dashed-dotted curve). The last case, which corresponds to completely neglecting electron-electron scattering in the simulation, yields two narrow gain spikes at around 2.8 and 3.6 THz, largely deviating from the experimental electroluminescence measurements.^{69} Ignoring the exchange effect (dotted curve) leads to an overestimation of the scattering, and thus a reduced peak gain and increased gain bandwidth. On the other hand, completely neglecting parallel spin collisions (dashed curve) leads to an underestimation of scattering and thus an enhanced peak gain.

### D. Elastic scattering processes

Impurity scattering has been shown to play an important role in QCLs,^{102,143,151} where it is often the dominant elastic scattering process. Also interface roughness can have a considerable impact on QCL operation.^{107,155–157} Alloy scattering occurs in semiconductor alloys such as AlGaAs and other ternary materials, and, usually, only has to be considered for the well material since the wave functions are largely localized in the wells.^{36} For elastic intrasubband scattering in the conduction band Γ valley, we have $k=k\u2032$ due to energy conservation, i.e., elastic intrasubband scattering does not affect the carrier distribution due to in-plane isotropy.

#### 1. Impurity scattering

Impurity scattering in QCLs arises from the doping, e.g., the ionized donors in the heterostructure. Charged impurities are to lowest order already considered in the Poisson equation, Eq. (20), corresponding to a mean field treatment of charges. An inclusion of higher order effects is in the semiclassical treatment obtained by corresponding scattering rates. The perturbation potential due to a charged impurity at position $(x\u2032,y\u2032,z\u2032)T$ is with $r\u2032=(x\u2032,y\u2032)T$ and $\u03f5=\u03f50\u03f5r,0$ given by the Coulomb potential

For the definition of the corresponding matrix element, various conventions exist.^{151,152} In analogy to Eq. (89) for electron-electron scattering, we define the bare matrix element as^{152}

The scattering rate is obtained from Eq. (57) as

with the form factor

Here, we have summed over the donors to include the effect of all ionized impurities, corresponding to an integration $S\u222b\u2026nD(z)dz$, with the doping concentration $nD(z)$. Furthermore, an additional factor *S*^{−2} has been added since the prefactor *S*^{−1} is omitted in Eq. (99). Due to energy conservation $Eik=Ejk\u2032$, the final wave vector magnitude $k\u2032$ is given by Eq. (59).

The total transition rate from a given initial state $|ik\u27e9$ to a subband *j* is found by summing over all final wave vectors $k\u2032$ using Eq. (61). We obtain^{54,152}

with $q(\theta )=|k\u2212k\u2032|=(k2+k\u20322\u22122kk\u2032cos\theta )1/2$. With Eq. (59), this yields

where $q02=2mj\u2225(Ei\u2212Ej)/\u210f2$. Values of θ resulting in complex $q(\theta )$ should be excluded from the integration in Eq. (102) since scattering is then impossible.

##### A. Screening

As for electron-electron scattering discussed in Sec. VI C, the random phase approximation can also be used to obtain the screened matrix elements $Vijs(z\u2032,q)$ for impurity scattering from the bare elements, $Vijb(z\u2032,q)$ in Eq. (99). We obtain^{151,152}

with $Vimjnb$ and $\Pi mn$ defined in Eqs. (89) and (94), respectively. In Eq. (100), $Vijb(z\u2032,q)$ has then to be replaced by $Vijs(z\u2032,q)$, and also Eq. (102) has to be changed correspondingly. As for electron-electron scattering, also here often simplified screening models using a constant screening wave number *q*_{s} are employed to avoid the numerical load associated with solving Eq. (104).^{54,151} Common models to handle screening are the model of Brooks and Herring^{158} and the model of Conwell and Weisskopf.^{159} While the model of Conwell and Weisskopf is based on a Rutherford-like scattering of electrons on bare Coulomb potentials with a cut-off radius, Brooks and Herring describe the impurity potential as being screened by free carriers. This latter approach is valid, if the constant screening length is much larger than the electronic wave length, which requires high temperatures and low carrier densities. For higher densities and low temperatures, the incorporation of more realistic screening has been shown to be essential,^{160} whereas the limit of negligible free carrier screening is better described by the approach of Conwell and Weisskopf.^{149} A detailed overview of the various refinements to the approach of Brooks and Herring has been given in Ref. 161. Electron scattering from screened impurities due to Brooks and Herring will be discussed in more detail in Sec. VIII D 3.

#### 2. Interface roughness scattering

It has been shown that scattering from rough interfaces can change the QCL performance significantly.^{107,155–157} There are two fundamentally different models for the interface roughness scattering in literature. In the first model, rough interfaces of quantum wells cause fluctuating well widths and fluctuating confinement energies of resonant states. In this model, the scattering potential is identified with the change of the well eigenenergies.^{162,163} In the second model, fluctuations of the band offset yields the scattering potential.^{164,165} The latter model may be favored for generality, since it does not require the existence of confined states, nor the distinguishability of level broadening by rough interfaces from other mechanisms. Nag has shown in Ref. 166 that both models agree well for quantum wells of various dimensions and materials. Thus, we consider interface roughness scattering due to the imperfections in the interface between the barrier and well material in the heterostructure, causing a local deviation of the interface $\Delta (x,y)=\Delta (r)$ from its average position *z*_{0} as illustrated in Fig. 14. Since it is not feasible to measure $\Delta (r)$ directly, the interface roughness is characterized by its standard deviation Δ and correlation length Λ. Typically, a Gaussian autocorrelation function is assumed^{162}

with $d=r\u2032\u2212r$. A change $\Delta (r)$ of the interface position corresponds to a perturbing potential

where *z*_{0} and *V*_{o} are the average interface position and the band offset, respectively. The “+” (“–”) sign is obtained if the barrier (well) is located at $z<z0$. Furthermore, *H* denotes the Heaviside step function. Since $\Delta (r)$ is small, i.e., on the order of a monolayer, we can approximate $\psi i,j(z0+\Delta (r))\u2248\psi i,j(z0)$. Thus, we obtain with Eq. (55)

and

with $q=k\u2212k\u2032$. Assuming a Gaussian autocorrelation Eq. (105), we obtain from Eq. (57) the scattering rate

where we additionally sum over all interfaces, located at positions *z _{n}*. Due to the energy conservation $Eik=Ejk\u2032$, the final wave vector magnitude $k\u2032$ is given by Eq. (59).

The total transition rate from a given initial state $|ik\u27e9$ to a subband *j* is found by summing over all final wave vectors $k\u2032$ using Eq. (61). We obtain

where $q(\theta )$ is again given by Eq. (103), and values of θ with $q2<0$ should be excluded from the integration since then scattering cannot occur.

The characterization of interface roughness by the parameters Λ and Δ is somewhat phenomenological, since it is difficult to directly extract them from experimental measurements of the interface profile. Rather, they are typically chosen to reproduce the measured optical transition linewidths.^{167} For GaAs-based terahertz QCLs, frequently used values are $\Lambda \u224810\u2009nm$ and $\Delta \u22480.1\u2009nm$.^{143,155} In Fig. 15, EMC simulation results for the gain profile are shown for two different types of terahertz QCLs.^{107} In addition to scattering of electrons with phonons, impurities, and other electrons, interface roughness scattering with different values of the interface roughness mean height has been considered. Results are shown for $\Delta =0$ (no interface roughness scattering), 0.12 nm and 0.3 nm (corresponding to one monolayer of GaAs). Depending on the assumed value of Δ, interface roughness can lead to considerable broadening of the gain profile and reduction in the peak gain. The phonon depopulation structure in Fig. 15(a) is less sensitive to interface roughness than the bound-to-continuum structure in Fig. 15(b), since the wave functions in the minibands of the latter structure are less localized and thus extend across more interfaces.

#### 3. Alloy scattering

In a ternary semiconductor alloy A_{x}B_{1–}_{x}C such as Al_{x}Ga_{1–}_{x}As, the A and B atoms are randomly distributed on the corresponding lattice sites. Thus, for a zinc-blende structure, each unit cell contains an atom of type A or B and an atom of type C. The potential in the semiconductor is then given by^{36,168}

where **x**_{n} denotes the position of unit cell *n*. Furthermore, *V _{n}* =

*V*

_{AC}(

*V*=

_{n}*V*

_{BC}) if the unit cell contains an atom of type A (B), where

*V*

_{AC}and

*V*

_{BC}are the potentials of the binary materials AC and BC. The conduction band profile of alloys is typically considered in virtual crystal approximation, obtained by averaging over the contributions of the two alloy materials AC and BC,

The perturbing potential corresponds then to the difference between the real conduction band energy $Vc(x)$ and *V*_{c,VCA} obtained by the virtual crystal approximation

with $\delta V(x)=VAC(x)\u2212VBC(x)$ and $cn=(1\u2212x)\u2009(cn=\u2212x)$ if the unit cell contains an atom of type A (B).

The scattering potential $\delta V$ is effective within the unit cell, and also screening effects are usually neglected due to the short-range nature of the scattering potential.^{36} Thus, on the scale of the envelope functions $\psi (x)$ which are slowly varying on the atomic scale, $\delta V$ can be modeled as a δ function potential^{36}

where $\Omega 0$ denotes the unit cell volume, given by $\Omega 0=a3/4$ for a zinc-blende structure with lattice constant *a*. Furthermore, $\delta V$ is the alloy scattering potential, which approximately corresponds to the conduction band offset of the binary materials involved, if they are not strongly lattice-mismatched.^{36} Generally, however, $\delta V$ deviates from the conduction band differences,^{156} and should be estimated from empirical tightbinding parameters.^{169–172} Using Eqs. (113) and (114), the matrix element Eq. (55) becomes with $q=k\u2212k\u2032$

and the square modulus is given by

For calculating the average of $|Vjk\u2032,ik|2$, the expectation value $\u27e8cncm\u27e9$ must be determined. Each unit cell has the probability *x* of containing an atom of type A, and 1 – *x* of containing an atom of type B. Thus, we obtain $\u27e8cncm\u27e9=0$ for $m\u2260n$ and $\u27e8cncm\u27e9=x(1\u2212x)$ for *m* = *n*. Furthermore replacing the sum over the unit cell positions $\u2211xnf(zn)$ by an integral $\Omega 0\u22121\u222bf(z)d3x=\Omega 0\u22121S\u222bf(z)dz$, we obtain

From Eq. (57), we obtain the scattering rate

where $\Omega 0(z),\u2009\delta V(z)$, and $x(z)$ have been taken into the integral to account for varying alloy compositions along the growth direction *z*. Due to the energy conservation $Eik=Ejk\u2032$, the final wave vector magnitude $k\u2032$ is given by Eq. (59). The total transition rate from a given initial state $|ik\u27e9$ to a subband *j* is with Eq. (61)

for $Ej<Eik$, and 0 otherwise.

### E. Self-consistent rate equation approach

In the self-consistent rate equation approach, the intersubband scattering rates in the rate equation Eq. (36) or Eq. (37) are self-consistently determined based on the corresponding Hamiltonian. Thus, this approach only relies on well known material parameters such as the effective mass, and not on experimental or empirical lifetimes as in Sec. V A. Offering a compromise between accuracy and predictive power on the one hand and relative numerical efficiency on the other hand, the self-consistent rate equation approach is widely used for the simulation of QCLs.^{67,77,91–93,173} Although extensions of the self-consistent rate equation approach to include the light field have been presented,^{77} the optical cavity field is neglected in most cases. Such simulations yield the unsaturated population inversion or gain, indicating if lasing can start at all, but giving no information about the actual lasing operation.

#### 1. Intersubband scattering rates

The number of electrons per unit energy and area in a thermalized subband *i* of a 2D system is given by $niEfiFD$ with the Fermi-Dirac distribution

where $EiF$ is here a “quasi” Fermi energy describing the kinetic energy distribution of the electrons within the subband *i*,^{132} and *T _{i}* is the associated electron temperature. Furthermore, $niE$ is the density of states per unit area and energy in a 2D system,

^{61}and $Eik=Ei+\u210f2k2/(2mi\u2225)$ is the energy of the electrons in subband

*i*with an in-plane wave vector

**k**. The scattering rate from an initial subband

*i*to a final subband

*j*is then obtained by averaging over the carrier distribution. Assuming a Fermi-Dirac distribution in the initial and final state, we obtain

^{132}

Here, $niE$ has been omitted in the enumerator and denominator because it is constant, $niE=mi\u2225/(\pi \u210f2)$ for $E\u2265Ei$. The index *m* denotes the corresponding scattering mechanism. $Wik,j(m)$ is, for example, given by Eqs. (82) and (88) for LO phonon and electron-electron scattering, respectively. The total rate is obtained by summing over the individual contributions of the different scattering mechanisms, i.e., $\tau ij\u22121=\u2211m(\tau ij(m))\u22121$.

The sheet density in a subband *i*, i.e., electron number per unit area, is given by

Furthermore expressing the carrier energy in terms of the in-plane wave vector, we obtain with $dE=\u210f2(mi\u2225)\u22121kdk$^{132}

An additional factor ${1\u2212fjFD[Ei+\u210f2k2/(2mi\u2225)\xb1E0]}$ can be included in the integral of Eq. (123), accounting for the final state blocking due to Pauli's exclusion principle.^{132} This correction is only relevant for high doping levels. Here, *E*_{0} is 0 for elastic scattering mechanisms and corresponds to the TO or LO phonon energy for optical phonon absorption (“+” sign) and emission (“–” sign), respectively. The Fermi energy can be calculated from Eq. (122)

For lightly doped semiconductors, as is often the case in QCLs, we have $(Eik\u2212EiF)\u226bkBTi$ in Eq. (120), which then approaches a classical Maxwell-Boltzmann distribution

Under this condition, Eq. (123) simplifies to

#### 2. Rate equations

In the self-consistent rate equation approach, typically the QCL is modeled as a biased periodic heterostructure, excluding effects such as domain formation.^{127} Then the simulation can be restricted to a single representative period far away from the contacts, additionally applying periodic boundary conditions.^{92} This corresponds to solving Eq. (36) with the self-consistently calculated scattering rates, Eq. (123). While Eq. (36) in principle includes the transitions to all equivalent levels in the different periods, in practice only scattering between the central period and adjacent periods has to be considered. Frequently the $112$ period model is used, which applies to QCLs where a period consists of an injector region and an active region.^{91,92} Here, for the active region states of the central period, only scattering transitions involving other states of the central period and the right-neighboring injector region are considered in Eq. (36). Analogously, for the states of the injector region, only scattering involving states of the central period and the left-neighboring active region are taken into account. Since the scattering rates Eq. (123) in general depend on the electron densities $nis$, a direct solution of Eq. (36) is not possible, and an iterative scheme is commonly used. For simulations without lasing included, i.e., $Wijopt=0$ in Eq. (36), setting $dtnis=0$ yields the steady state solution^{92}

where $i=1..N$ refers to the central period containing *N* subbands. The summation index *j* now only includes subbands in the central period and adjacent periods, as discussed above. For subbands outside the central period, the sheet density $njs$ of the equivalent level in the central subband is used. For the numerical computation of the subband populations, initially identical electron densities $nis=ns/N$ are assumed. Each iteration involves computing the scattering rates with Eq. (123) or Eq. (126), calculating the new values for $nis$ using Eq. (127), and renormalizing the sheet densities so that Eq. (38) is fulfilled. Convergence can be accelerated by combining the sheet densities of the previous two iterations, $\xi nis,new+(1\u2212\xi )nis,old$, as input for the next iteration,^{91} where the relaxation parameter is typically chosen as $\xi =0.5$. The simulation has converged when the obtained sheet densities do not significantly change anymore between iterations.

#### 3. Kinetic energy balance method

Since *T _{i}* in Eq. (120) is unknown, the scattering rates [Eq. (123)] have to be evaluated assuming that the electron temperature in each subband is equivalent to the lattice temperature,

*T*=

_{i}*T*

_{L}. However, the electron temperature can significantly exceed

*T*

_{L}in quantum cascade lasers.

^{174}Thus, the rate equation model has been extended to account for electron heating, typically assuming an identical electron temperature

*T*=

_{i}*T*

_{e}for all subbands.

^{173}The kinetic electron energy generation rate per period and unit device in-plane cross section is given by

^{173}

where *i* sums over the subbands of the central period and *j* also includes states from adjacent periods, as discussed above. Furthermore, *m* sums over the different scattering contributions, such as LO phonon emission $(E0(m)=\u2212ELO)$ and absorption $(E0(m)=ELO)$, elastic scattering mechanisms $(E0(m)=0)$ and electron-electron scattering. For inelastic scattering mechanisms such as phonon-induced transitions, also intrasubband contributions *i* = *j,* have to be included.^{173} Electron-electron scattering is often neglected in Eq. (128) since the net kinetic energy does not change for the important cases where both electrons stay within their respective subbands or swap the subbands. Also photon transitions do not change the kinetic energy for parabolic subbands due to **k** conservation.^{77} The numerical evaluation proceeds as follows. The average electron temperature corresponds to the value of *T*_{e} where *R ^{E}* = 0 in Eq. (128). The rate equation is self-consistently solved as described above for an initial guess of

*T*

_{e}, e.g.,

*T*

_{e}=

*T*

_{L}. Based on the obtained $nis$ and $\tau ij(m)$,

*R*in Eq. (128) is calculated. This procedure is repeated and the guess for

^{E}*T*

_{e}is iteratively improved until $RE\u22480$ is obtained. Since the assumption of a single effective electron temperature is not always adequate in QCLs,

^{175}extended approaches have been developed allowing for different effective temperatures

*T*in the individual subbands.

_{i}^{93}

### F. Ensemble Monte Carlo method

Semiclassically, the carrier transport between states $|ik\u27e9$ in a quantum well system is given by the Boltzmann equation^{176}

with the scattering rates $Wik,jk\u2032=\u2211mWik,jk\u2032(m)$, where *m* sums over the individual contributions of the different scattering mechanisms. The distribution function $fik(t)$ represents the probability of the state $|ik\u27e9$ being occupied at a given time *t*. Eq. (129) corresponds to an extended version of the rate equations Eq. (36), where *i* and *j* have been replaced by *i***k** and $jk\u2032$, respectively. The physical quantities of interest, such as the sheet densities $nis$ and current density can be extracted from $fik(t)$. For the numerical evaluation of Eq. (129), typically the EMC method is used, which is based on statistical sampling of the scattering events for a large ensemble of carriers,^{149} here $Ne\u2248104..105$ electrons. The large number of carriers considered allows for an extraction of the physical quantities as function of *t* by statistical averaging over the carrier ensemble. Thus, this method is also applicable to time dependent processes where we cannot use temporal averaging as for ergodic systems.^{149} Furthermore, electron-electron scattering can be implemented as a two-electron process where a second electron is randomly chosen as scattering partner.^{141}

#### 1. Simulation technique

Figure 16 contains a schematic diagram of the EMC algorithm. The system dynamics is evaluated up to a time *t*_{sim}, which must be chosen long enough to ensure convergence to the stationary solution. The simulation is divided into subintervals $\Delta t$, where the scattering dynamics is subsequently evaluated. $\Delta t$ should be chosen so small that the average electron distribution does not change significantly over the time interval, but big enough to include several scattering events per electron on average. Assuming periodic boundary conditions, we can restrict our simulation to a few periods. Each electron is characterized by its subband *i* and in-plane wave vector **k**. All rates $Wik,j(m)$ for the various scattering mechanisms *m* are computed and tabulated at the beginning of the simulation to save computational resources. Here, we have to introduce a discrete grid for the wave vector. Since in the conduction band Γ valley, all states $|ik\u27e9$ in subband *i* with the same value $k=|k|$ are equivalent due to in-plane isotropy, it is practical to use the kinetic energy $Ekin=\u210f2k2/(2mi\u2225)$ instead. The kinetic energy grid then divides the energy axis into segments *n* of widths $\Delta E(n)$ centered around discrete energies $Ekin(n)$ with $k(n)=(2mi\u2225Ekin(n))1/2/\u210f$. However, some rates $Wik,j(m)$, such as Eq. (88), depend on the initially unknown carrier distribution itself, i.e., $Wik,j(m)=Wik,j(m)(t)$. This problem can be overcome by tabulating an upper estimate $W\u0303ik,j(m)$ for time dependent scattering rates and compensating for the too high value by introducing artificial “self-scattering,” as described further below.

An important quantity is the carrier distribution function, which can directly be obtained from $Ni(n)$, denoting the number of simulated electrons in the *n*th energy cell of subband *i* at a time *t*. With the density of states per unit area and energy in a 2D system $niE=mi\u2225/(\pi \u210f2)$,^{61} the number of available states in the *n*th energy cell is $niE\Delta E(n)S$. The simulated device in-plane cross section is with the sheet doping density per period *n*^{s} given by $S=Ne/(nsNp)$, where *N*_{p} corresponds to the number of periods over which the *N*_{e} simulated electrons are distributed. The carrier distribution function is then approximately given as

where *n* indicates the energy cell containing the value *E*_{kin}. Equation (130) can also be expressed as a function of $k=(2mi\u2225Ekin)1/2/\u210f$. Additional temporal averaging of the carrier distribution function Eq. (130), for example over a simulation subinterval $\Delta t$, further reduces the fluctuations resulting from the treatment of the carriers as discrete particles.

The simulation is initialized by assigning values *i* and **k** to each electron in the ensemble. Here, in principle arbitrarily chosen values can be used since Eq. (129) will converge to its stationary solution after a sufficiently long simulation time *t*_{sim} independently of the chosen initial conditions. However, the required *t*_{sim} to obtain convergence can be reduced by assuming a suitable initial electron distribution which is not too far from the converged solution. For example, the carriers can be distributed equally within the subbands, and a thermalized distribution Eq. (125) in each subband can be assumed. The Monte Carlo method is based on a stochastic evaluation of the scattering events, assuming that the electron undergoes scattering after a randomly selected free flight time *t*_{f}. For a time independent scattering rate $\tau 0\u22121$, the free flight time is given by $tf=\u2212\tau 0ln(r)$, where *r* is a random number evenly distributed between 0 and 1.^{149} Here, we choose $\tau 0\u22121\u2265\u2211m\u2211jW\u0303ik,j(m)$, corresponding to an upper bound for the total outscattering rate $Wik(t)=\u2211m\u2211jWik,j(m)(t)$, where *j* = *i* is included in the summations to account for intrasubband scattering. For each ensemble electron, the scattering dynamics is evaluated over a simulation subinterval $\Delta t$, where the last free flight is continued in the next subinterval. After each scattering event, the final subband *j* of the scattered electron and scattering mechanism *m* are randomly selected based on their associated probability $\tau 0Wik,j(m)(t)$. The too high amount of scattering assumed is corrected by introducing artificial “self-scattering,” which occurs with a probability $1\u2212\tau 0\u2211m\u2211jWik,j(m)(t)$ and does not change the carrier state at all.^{149} Pauli's exclusion principle can be considered by subsequently rejecting the occurred scattering event with a probability $fjk\u2032$. The final wave vector magnitude $k\u2032$ is obtained from Eq. (59) for elastic scattering, and from Eq. (60) for inelastic processes. Electron-electron scattering is in EMC typically implemented as a two-electron process with a randomly selected partner electron and scattering angle. To conserve energy for each scattering event individually, also the partner electron has to undergo scattering, i.e., the final wave vector has to be determined for both participating electrons.^{141}

The periodic boundary conditions can, for example, be implemented by simulating three periods, where the *N*_{e} electrons are located in the central period, i.e., *N*_{p} = 1 in Eq. (130). Electrons scattered to the first or third period are automatically injected into the equivalent subband of the central period.^{99} By counting the difference $\Delta N$ of electrons scattered from the central period to the left- and right-neighbouring period, respectively, over a simulation subinterval $\Delta t$, the current density can be computed as

### G. Inclusion of the optical cavity field

The optical gain and photon-induced transition rate are given by Eqs. (52) and (54) for a transition from a level *i* to a level *j*, which here correspond to states $|ik\u27e9$ and $|jk\u2032\u27e9$, respectively. The dipole matrix element then becomes

with

and for $|dik,jk\u2032|\u20092$, we obtain

Here, we have used that $|\u222bSexp[i(k\u2032\u2212k)r]d2r|2$ can be approximated by $4\pi 2S\delta (k\u2032\u2212k)$ for sufficiently large in-plane cross sections *S*. The photon-induced transition rate given in Eq. (54) now becomes

where we sum over all relevant cavity modes with frequencies $\omega m$ and intensities *I _{m}* to account for multimode lasing. Adapting Eq. (51) to the present case, the definition of the Lorentzian lineshape function becomes

where $\omega ij(k)=(Eik\u2212Ejk)/\u210f$ denotes the resonance frequency and $\gamma ij(k)$ is the optical linewidth of the transition, which is given by

when only the lifetime broadening contributions are considered.^{107} The total photon-induced transition rate from a given initial state $|ik\u27e9$ to a subband *j* is found from Eq. (135) by summation over all final wave vectors $k\u2032$ using Eq. (61)

The transition rate due to spontaneous photon emission can also be directly calculated from Eq. (138),^{177} but is usually negligible compared to other scattering mechanisms in QCLs. The gain contribution at frequency ω of a single electron in state $|ik\u27e9$, i.e., $nis=1/S$, is obtained from Eq. (52) by summing over the transitions to all available final states $|jk\u2032\u27e9$. With Eqs. (61) and (134), we obtain

For $Eik<Ejk$, we have $\omega ij(k)/|\omega ij(k)|=\u22121$, indicating absorption and thus resulting in a negative contribution to the optical gain.^{76} For EMC simulations, the gain can be evaluated by summing Eq. (139) over all simulated electrons *n* in the corresponding states $|ink\u27e9$ of the central simulation period(s)^{76}

where *j* also includes states in adjacent periods. *S* then corresponds to the simulated device in-plane cross section, which is for a sheet doping density per period *n*^{s} given by $S=Ne/(nsNp)$, where *N*_{P} corresponds to the number of periods over which the $Ne$ simulated electrons are distributed.

For rate equations, the electron wave vector dependence is not considered, and thus a *k* independent averaged value $\gamma ij$ has to be taken in Eq. (136). Furthermore, $\omega ij$ does not depend on the wave vector if nonparabolicity effects are neglected. Then the gain simplifies to Eq. (52), thus becoming

The intensity evolution for a mode *m* is again given by Eq. (49)

The carrier transport and intensity evolution have to be simulated by a coupled approach, as illustrated in Fig. 17. The carrier transport simulation based on rate equations or EMC now includes the photon-induced transition rates Eq. (138), which depend on the intensities given by Eq. (142). On the other hand, Eq. (142) depends on the gain which is extracted from the carrier transport simulations by using Eq. (140) or Eq. (141).

### H. Selected EMC simulation results

In the following, some EMC simulation results are presented for both mid-infrared and terahertz QCLs, illustrating the versatility of this approach. Examples without inclusion of the optical cavity field are shown in Secs. VI A–VI D, see Figs. 11, 13, and 15 where results for the unsaturated optical gain are displayed. Here, we focus on coupled simulations of the carrier transport and the optical cavity field.

The EMC method has been applied to study the optical power and gain saturation in terahertz QCLs, yielding good agreement with experimental data.^{76} Here, examples are presented for a high temperature QCL lasing up to 164 K.^{178} In Fig. 18, the simulated unsaturated spectral gain curve is compared to the result obtained for carrier-light coupling included, demonstrating gain saturation and clamping at the threshold gain value due to the lasing field. In Fig. 19, the obtained optical output power as a function of lattice temperature is compared to experimental data, demonstrating good agreement with experiment.

Due to the stimulated optical transitions, the lasing field significantly influences the carrier transport, affecting not only the subband populations and thus the optical gain but also the electric current. This is especially the case for high efficiency mid-infrared QCLs. In Fig. 20, the measured current density of a mid-infrared QCL with ≈50*%* wall-plug efficiency^{9} is compared to EMC results with and without carrier-light coupling included.^{106} Good agreement with experiment is only obtained in the former case, while the simulation neglecting photon-induced processes significantly underestimates the current density, demonstrating the significant influence of these processes on the carrier transport.

## VII. DENSITY MATRIX APPROACHES

The one-dimensional density matrix approach is frequently used for the analysis and optimization of QCLs.^{17,111–113} Also three-dimensional versions have been developed.^{114,115} One-dimensional density matrix approaches can be seen as a quantum mechanical generalization of rate equations, Eq. (36), to include effects such as resonant tunneling and dephasing. In the semiclassical description, transport through a barrier occurs instantaneously due to electrons being scattered into wave functions which are spatially extended across the barrier,^{125} see Fig. 21(a). In the quantum mechanical picture, the electron transport across the barrier is described by a coherent superposition of the extended states with a narrow anticrossing energy gap $\Delta E$, resulting in a localized electron wavepacket. Due to the coherent time evolution of these states, the wavepacket oscillates between the left and right well with the so-called Rabi oscillation frequency $\Omega =\Delta E/\u210f$,^{125} and the corresponding tunneling time corresponds to half the duration of an oscillation cycle, $\tau tun=\pi /\Omega $.^{111} The semiclassical picture is adequate as long as the Rabi oscillations are not significantly dampened by dephasing, since then the wavepacket oscillates uniformly between the wells and the averaged population distribution corresponds to the semiclassical description. However, for strong dephasing with relaxation times $\tau \u2264\tau tun$, the wave packet no longer oscillates. Rather, the electrons accumulate in the left well, and the semiclassical picture of a uniform electron distribution across the extended states is no longer valid.^{125} Now the current transport is not limited by the scattering-induced electron transport, but by tunneling through the barrier. This effect is commonly referred to as wave function localization due to dephasing, and occurs for thick barriers with narrow anticrossing gaps $\Delta E$ and thus long tunneling times $\tau tun$.^{125} Especially in terahertz QCLs, thick injection barriers are used with $\Delta E\u22481meV$,^{125} corresponding to $\tau tun$ of around 2 ps, while the dephasing time can be estimated from the spontaneous emission linewidth of ∼5 meV (Ref. 125) to be ≈0.3 ps.

Rather than describing localized wavepackets by a coherent superposition of eigenstates, frequently localized wave functions are used,^{111–113} see Fig. 21(b). Here, the wave functions are computed for each QCL period separately, by assuming that the thick injection barriers at the left and right sides of the period are infinitely thick.^{111} The Schrödinger equation, Eq. (4), is then solved for the corresponding tight-binding conduction band profile *V*_{tb} rather than the actual potential *V*. The corresponding Hamiltonian is obtained in the framework of tight-binding theory.^{36} Specifically, the Rabi frequency for a doublet of states spanning the coupling barrier is given by $\Omega ij=\u210f\u22121\u27e8i|V\u2212Vtb|j\u27e9$.^{111} Depending on the investigated QCL structure, resonant tunneling may have to be considered at more than a single coupling barrier in each period.^{112,113} Then the period is subdivided into smaller regions, and the tight binding formalism has to be applied accordingly.

Wave functions $\psi i$ represent pure states, and thus are not adequate for describing decoherence effects due to the interaction with the environment. This can be accomplished by means of a density matrix $\rho ij$ (see also Sec. V B), which can represent both pure and mixed states and thus is suitable to model dephasing. The diagonal terms $\rho ii$ correspond to the occupation of level *i*, and the off-diagonal elements $\rho ij$ are the coherence or polarization terms for the doublet *i, j*.^{125} The density matrix can be normalized so that the diagonal terms correspond to the sheet density, $\rho ii=nis$. The time evolution of the density matrix $\rho ij(t)$ is described by the von Neumann equation

where the Hamiltonian matrix elements are defined as $Hij=\u27e8i|H\u0302|j\u27e9$. For describing the carrier transport in QCLs, the Hamiltonian is divided into two parts, $H\u0302=H\u03020+H\u0302\u2032$. Here, $H\u03020$ describes the coherent evolution of the quantum system due to the conduction band potential, while $H\u0302\u2032$ contains the perturbation potentials corresponding to the various scattering mechanisms and introduces dissipation to the system. $H\u0302\u2032$ is commonly implemented in a somewhat phenomenological manner by using transition and dephasing rates, similar as for the rate equation approach discussed in Secs. V A and VI E. Equation (143) can then be cast into the form^{111,113}

where $\omega ij=(Ei\u2212Ej)/\u210f$, and $\tau i\u22121=\u2211j\u2260i\tau ij\u22121$ indicates the total inverse lifetime of level *i*. The scattering rates $\tau ij\u22121$ have already been discussed in Secs. V A and VI E. The dephasing rate is given by^{125}

Here, $(\tau i\u22121+\tau j\u22121)/2$ is the lifetime broadening contribution [see also Eq. (137)], and $\tau pure,ij\u22121$ contains the pure dephasing. The Rabi frequency $\Omega ij$ is nonzero only for doublets spanning a coupling barrier; specifically, $\Omega ii=0$. In Eq. (144), resonant tunneling is assumed to be independent of the in-plane wave vector, and the rates $\tau ij\u22121,\gamma ij$ are averaged over the kinetic electron distribution within the subbands.^{111} The stationary solution of Eq. (144) is obtained by setting d_{t} = 0, yielding^{113}

with

Furthermore, Eq. (38) has to be fulfilled, i.e., the total sheet density in each period is determined by the doping sheet density *n*^{s}.

Equation (144) can be solved using empirical rates, or implemented in a self-consistent manner by calculating the scattering rates $\tau ij\u22121$ as discussed in Sec. VI E. The pure dephasing rate $\tau pure,ij\u22121$ in Eq. (145) can be computed based on intrasubband scattering transitions,^{113,179} but is often treated by assuming an empirical value.^{112,125} In a simplified model, the transport across the barrier can be restricted to the coupling between two states 1 and 2. For this case, Eq. (146) simplifies with $\tau 12\u22121=\tau 21\u22121=0$ to $\tau 2\u22121n2s+R12(n2s\u2212n1s)=0$. The current density through the barrier is then given by $J=eR12(n1s\u2212n2s)=e\tau 2\u22121n2s$. With Eq. (147) and $n1s+n2s=ns$, we obtain^{10}

In the three-dimensional density matrix method, additionally the in-plane wave vector **k** is taken into account, i.e., the electrons are described by states $|ik\u27e9$ rather than $|i\u27e9$. Consequently, the density matrix is given by $\rho ik,jk\u2032$, and $i,j,\u2113$ in Eq. (143) have to be replaced by $ik,\u2009jk\u2032$, and $\u2113k\u2032\u2032$, respectively. Three-dimensional density matrix approaches can be seen as a quantum mechanical generalization of the Boltzmann equation given in Eq. (129). Various approaches based on the three-dimensional density matrix have been developed for QCL simulation, where the scattering mechanisms are self-consistently implemented based on the corresponding Hamiltonians.^{114,115} Specifically, also a hybrid density matrix-Monte Carlo approach has been introduced, where the tunneling transport through the coupling barrier is treated based on the density matrix formalism, while scattering inside each period is treated semiclassically using an EMC approach.^{125}

## VIII. NON-EQUILIBRIUM GREEN'S FUNCTION FORMALISM

### A. General scope of the non-equilibrium Green's function method

Quantum cascade devices utilize charge transport in structures of the nanometer length scale. In this regime, quantum effects such as coherent tunneling, interference, and confinement play a very important role for the transport physics. At the same time, however, the devices are run at finite temperatures, which support a significant amount of scattering with phonons. Since most devices are doped, alloyed and/or based on heterojunctions of different materials, impurity and alloy disorder as well as surface and interface roughness influence the transport, too. All these scattering effects share the fact that the full quantum information (i.e., the full phonon phase, the precise position of the impurities, etc.) is either unknown or lost in the statistics of large numbers of scattering events. Consequently, these effects contribute to incoherent scattering and dephasing.

It is well established that the NEGF theory is the most general scheme for the prediction of coherent and incoherent quantum transport. Since its introduction in the 1960s,^{180–182} this formalism has been successfully applied on a great variety of transport problems. These problems include, but are not limited to spin,^{183–185} phonon,^{186–188} and electron transport,^{189,190} covering different materials such as metals,^{191,192} semiconductors,^{193} topological insulators,^{194} and even various dimensionalities such as layered structures,^{195} nanotubes,^{187,196} fullerenes,^{197} and molecules.^{198,199} For electronic transport, the NEGF method has been implemented in a large variety of representations, ranging from envelope function approximations such as the effective mass^{200,201} and **k·p** method,^{202} to atomistic representations such as tight binding^{203,204} and even density functional theory models.^{205} In most cases, the NEGF method is used for transport in open systems, i.e., devices connected to spin, charge or heat reservoirs via semi-infinite leads. Nevertheless, two different approaches have been successfully applied to mimic the field-periodic conditions cascade structures are facing.^{118,122,206,207}

It is this high flexibility of NEGF that offers choosing specialized basis functions such as, e.g., Wannier-Stark states for cascade structures,^{206} wire modes for homogeneous wire systems,^{189} or generic basis functions for general cases within the low rank approximation.^{208} It is true for all these cases that the closer the basis functions match the actual quasi particles of the system, the fewer basis functions are required to reliably predict the device performance and the more efficient the numerical implementation of NEGF will be. It is worth to mention here that the NEGF method natively contains more information than semiclassical methods such as the Boltzmann equation. In fact, it has already been shown in the 1960s how to approximate the NEGF equations to yield the Boltzmann equation.^{181}

The NEGF method, however, faces one major drawback compared to most other methods: the numerical solution of the NEGF equations is expensive both in terms of memory and central processing unit (CPU) time. The numerical costs are particularly high if incoherent scattering in the self-consistent Born approximation is considered. For this reason, most numerical implementations of NEGF on real devices do not include electron-electron scattering beyond the first order (Hartree) approximation. Higher orders of electron-electron scattering (i.e., exchange terms) are required to model energy transfer during inelastic electron-electron scattering, but their numerical load typically prohibits the implementation on concrete transport problems.

### B. Overview of the non-equilibrium Green's function method

#### 1. Fundamental equations and observables

As discussed in Sec. II A, electrons in quantum cascade lasers can be successfully described within the effective mass approximation. Hereby, the devices are typically considered as laterally homogeneous quantum well heterostructures. The electronic structure is then represented with the Hamiltonian

where *k* is the lateral electron momentum, $\Phi (z)$ the electrostatic potential, and $Vc(z)$ denotes the material and position dependent conduction band edge, including the band offsets. Note that the effective mass is energy-dependent to include the nonparabolicity as described in detail in Sec. II C. In the stationary limit, the NEGF method describes transport with the electronic retarded and lesser Green's function $GR$ and $G<$, respectively. These functions solve four coupled partial differential equations that read in operator form for a given energy *E* and in-plane momentum **k**^{209}

*D*and

^{R}*D*

^{<}are the sum of retarded and lesser Green's function of the environment. Equations (150a) and (150b) are also referred to as Dyson and Keldysh equation, respectively. The expressions for the self-energies, Eqs. (150c) and (150d), are discussed in Sec. VIII D for various scattering mechanisms. We note that in that section and in the rest of this paper all Green's functions and self-energies are given in real space representation which requires a transformation of Eqs. (150) from operator space into real space. Thus, the Green's functions and self-energies given in the following have the dimensions eV

^{−1}m

^{−1}and eV m

^{−1}, respectively, unless stated otherwise. Green's functions and self-energies are solved in the self-consistent Born approximation, i.e., they are solved iteratively until convergence is achieved. Most of the known current conserving simplifications of the self-energies do not ease the convergence of the self-consistent Born iterations, but reduce the numerical burden of each individual iteration only. In contrast to methods that require the solution of the Schrödinger equation, the solutions of Eqs. (150) do not require to solve an eigenvalue problem. Consequently, energy dependent effective masses in the Hamiltonian $H\u03020$ do not increase the numerical complexity of NEGF. The Green's functions and self-energies are functions of two spatial coordinates $z,z\u2032$, the absolute value of the lateral momentum

*k*and the energy

*E*. The energy and spatially resolved spectral function $A(z,E)=i[GR(z,z,0,E)\u2212GR\u2020(z,z,0,E)]$ shows width and location of resonant states in the system. This indicates that

*G*contains the information of resonant states and the density of states. The

^{R}*G*

^{<}function contains in addition information of the occupancy of the electronic states. Consequently, occupancy related observables such as density, current, and optical gain are dependent on

*G*

^{<}: the spatially and energy-resolved density $n(z,E)$ and current density $J(z,E)$ are defined in relation to the density $n(z)$ and current density $J(z)$, respectively,

The current density in Eq. (152) is only correct as long as the kinetic energy operator is the only term in the Hamiltonian that does not commute with the position operator. This can be seen from the fundamental probability current density operator in real space **x** and time *t* representation (see, e.g., Refs. 210 and 211)

with the velocity operator

Please note that the Green's function in Eq. (153) is given in real space **x** and time *t* representation and has the dimension eV^{−1} m^{−3} s^{−1}. The optical field amplitude absorption coefficient $\alpha (z,\omega )$ for a photon of frequency ω is a function of the permittivity $\u03f5(z,\omega )$^{212}

Here, $\u03f50$ is the vaccum permittivity and *c* denotes the speed of light. The optical absorption coefficient is usually defined with respect to the field intensity (*a*) rather than with respect to the field amplitude,^{88} see Sec. III. In particular, in the gain regime one often refers to the power gain −*a*. These quantities are related by $a=2\alpha $. The permittivity depends on the complex conductance $\sigma (z,\omega )$ and the materials dielectric constant $\u03f5r(z)$

Before lasing starts and the perturbation $\delta V\u0302(\omega )$ due to the optical field is still small, the optical absorption can be extracted from the linear response of the Green's functions. To first order, the change of the lesser Green's function is given by (for a given **k**)^{116}

Changes to the scattering self-energies $\delta \Sigma $ due to the photonic perturbation $\delta V\u0302$ are first order vertex corrections.

The wavelengths of terahertz lasers are much larger than the typical dimensions of the QCL periods. Thus, for the electron transport calculation, the optical electric field can be considered constant in the active device. Since the QCL laser light is usually linearly polarized in transport direction *z*, the perturbing potential reads in Coulomb gauge

with the photon electric field component in *z* direction *E _{z}*. This results in the change of the current density $\delta J(z,\omega )$ in linear order of $\delta V\u0302(\omega )$

The first term of the last equation results from the current operator in Eq. (152) applied on the change of the lesser Green's functions $\delta G<$. The second term results from the change of the velocity and of the current operator due to the fact that $\delta V\u0302(\omega )$ does not commute with the position operator [see Eq. (154)]. The quotient of the perturbation of the current density $\delta J(z,\omega )$ and the electric field $Ez(\omega )$ of the photon gives us the optical conductance

If vertex corrections $\delta \Sigma $ are ignored for simplicity, Eqs. (158), (157), and (159) can be combined into an equation for the optical conductance

When the self-consistently solved Green's functions are used for the optical conductance, Eq. (161) fully accounts for the self-consistently calculated electron states and their non-equilibrium state occupations. The vertex corrections $\delta \Sigma $ in Eq. (157) increase the numerical load significantly, since they require self-consistent iterations of $\delta G$ with $\delta \Sigma $. Both $\delta G$ and $\delta \Sigma $ depend on the optical frequency ω in addition to the “standard” dependence of the Green's function [e.g., $\delta \Sigma =\delta \Sigma (z,z\u2032,k,E,\omega )$]. Thus, such self-consistency is expensive in terms of memory and time. It has been shown that vertex corrections narrow the optical linewidths and increase the peak height of the absorption coefficient in terahertz QCLs.^{117} Calculations without $\delta \Sigma $ may not predict quantitative values of $\alpha (z,\omega )$, but can serve for qualitative predictions only. If the interaction with the photonic field should be considered beyond linear response—such as in the case of electron transport during lasing, time dependent NEGF with a periodically oscillating electric field is required. Details of this approach can be found in Ref. 124.

The balance between the benefits of the NEGF method and its numerical load is strongly device physics dependent. If the electron transport is clearly dominated by incoherent scattering and tunneling across multiple barriers is negligible, semiclassical models are clearly more efficient than NEGF. Otherwise, if scattering is negligible, the NEGF method has to compete with numerically more efficient methods such as the Schrödinger equation or the quantum transmitting boundary method.^{213} If scattering in low order captures the physics and energy resolved information is not desired, the density matrix method is more efficient than NEGF as well. The strength of NEGF is that it allows to predict transport in all these regimes and gives deep insight to the ongoing processes in any of the before mentioned situations. It will be discussed in Sec. VIII E in detail that a single terahertz QCL can move from the ballistic to the scattering dominated regime with the applied electric field. Therefore, it is appropriate to use NEGF on these devices.

#### 2. Different basis representations and low rank approximations

The solution of the NEGF equations, Eqs. (150), requires many matrix-matrix products and matrix inversions. Both scale cubically with *N*, the number of degrees of freedom the equations are discretized in. As a result, highly resolved NEGF calculations can easily require modern supercomputers.^{44} It is obvious that any reduction of *N* will reduce the number of floating point operations signficantly.^{208} Real space discretization of Eqs. (150) faces the challenge that any device feature such as differing widths of semiconductor layers has to match to the chosen resolution. Particularly small device features require an inhomogeneous grid in real space, increasing the numerical complexity.

The most common technique to reduce *N* is to transform the NEGF equations into a system of basis functions that are equal or close to the actual propagating eigenstates of the system's Hamiltonian. The more the basis functions agree with the eigenfunctions of the system, the less the basis functions couple with each other and the fewer functions need to be included in the actual transport calculation. If the target basis has less states than the rank *N* of the original (real space) representation, the transformation matrices are rectangular. In the ideal case, when the basis functions are isolated from each other, their individual contribution to transport can be solved and summed. Prominent examples of such rank reductions are the mode space approaches for transport in nanowire structures. If the basis functions are coupled, it is common to still limit the number of considered functions. In this case, however, this approach gives only approximate results. In general cases, transformations with rectangular matrices that reduce the rank of the equations are called low rank approximations.^{208}

In the framework of cascade devices, a small number *N* of Wannier or Wannier-Stark basis functions have proven to represent quantum transport well enough for reliable predictions. For instance, Lee *et al.* were able to solve the IV characteristics and optical output performance of terahertz QCLs with only 5 states per period.^{214} Although it is a common and numerically very efficient technique, caution is required when results of the low rank space are transformed back into the original real space representation if the rank of the two differs a lot. Numerically, the back transformation produces results with the high resolution of the original space, but the reliability of the data in the real space resolution strongly depends on how well the basis functions of the low rank space match the physics of the device. Unphysical oscillations of the current density with real space can indicate such reliability issues. Increasing the number of considered basis functions in the low rank space can improve the reliability of the approximate results.^{208,214}

### C. Boundary conditions for cascade systems

Boundary conditions in stationary NEGF always fall into two categories: boundary conditions for *G ^{R}* that are related to the density of states and boundary conditions for

*G*

^{<}that are related to state occupancy and particle distributions. Since quantum cascade lasers consist of many repetitions of the same building block—a sequence of quantum wells and barriers—NEGF is often applied on cascade structures with field periodic boundary conditions.

^{121,206,215}Field periodic boundary conditions for the retarded Green's function can be mimicked by solving the NEGF equations in a system of 3 or more periods: the innermost period is considered as the actual active device while the boundary periods are considered as charge suppliers only. The Poisson equation is solved with periodic boundary condition and an applied homogeneous electric field is added to the resulting potential. Since the source sided periods are on average on a higher potential (and contrary for the drain sided periods), the system effectively faces field periodic charge boundary conditions. The initial guess for the lesser Green's function

*G*

^{<}is also field periodic. Since the system is closed, self-consistent Born iterations of the higher order scattering terms may change the number of electrons. To maintain global charge neutrality, the

*G*

^{<}function needs rescaling between iterations.

^{215}Although this type of boundary treatment is appealing for cascade systems, care on the numerical details has to be taken: In the first iteration of the self-consistent Born approximation the system is purely coherent. Then, the equation for

*G*is van-Hove singular whenever the energy

^{R}*E*agrees with a resonance of the system. Otherwise, the resulting retarded Green's function is real and the spectral function

*A*vanishes (i.e., no electronic states are found). To solve this issue, a small artificial retarded scattering self-energy is typically included in Eqs. (150). This self-energy maintains a minimal broadening of resonances. This artificial broadening has to be negligible compared to the actual scattering self-energies. Note that purely coherent and field periodic calculations would clearly yield divergent electron energies in nonequilibrium.

Most applications of NEGF tackle transport in open systems. In open systems, the active device is connected via semi-infinite leads with charge reservoirs that supply or collect electrons with a constant, typically equilibrium distribution. Technically, the boundary conditions of *G ^{R}* and

*G*

^{<}are then introduced by including contact self-energies $\Sigma conR$ and $\Sigma con<$ in the equations for

*G*and

^{R}*G*

^{<}in Eqs. (150), respectively. In this case, the large number of cascade periods can be mimicked by considering one or a few periods as the actual active device and including the potential landscape and material sequences of the adjacent periods in the calculation of $\Sigma R$.

^{118,122}In this way, the presence of the leads broadens resonant states in the active device and allows for the solution of purely ballistic transport without adding artificial retarded self-energies. It also allows for the assessment of individual scattering mechanisms, since scattering and its broadening of the states is not essential for numerical convergence and can be turned off selectively. However, the boundary condition for

*G*

^{<}requires a model for the electron distribution in the leads. It is common in NEGF for open systems to assume equilibrium Fermi distributions. In the case of cascade devices, however, such equilibrium distributions do not include possible heating effects within previous cascade periods and therefore allow results of nonperiodic electron distributions: When heating effects exist, the electrons do not dissipate all energy they gain while traversing the potential drop of a single period. Then, the electron distributions at source and drain side of each period differ.

^{118}

Alternatively, an electron distribution can be read out of *G*^{<} at device positions that are equivalent to the lead/device interfaces. In this case, the nonequilibrium distribution is defined as

with the spectral function

In equilibrium, the function $f(z,k,E)$ equals the Fermi distribution $f(E,\mu )$.^{209} When the distribution function of source (*z* = 0) and drain (*z* = *L*) are assumed to be equivalent to the distribution functions at positions in the device that are a single QCL period apart

field periodic boundary conditions for *G*^{<} are mimicked.^{122} Here, *L*_{p} is the length of a QCL period, $\varphi $ is the potential drop per period, and *n* is the number of explicitly considered QCL periods. To ensure global charge neutrality, the Poisson equation is solved with the applied electric field at both boundaries.

### D. Scattering self-energies

In the following, the resulting expressions for the self-energies, given in Eqs. (150c) and (150d), are summarized for the most important scattering mechanisms in QCLs. If multiple scattering mechanisms are considered, the individual self-energies of the respective scattering mechanisms can be summed up and the summed self-energy is then used in Eqs. (150). Some of the scattering self-energies are numerically so expensive to implement that approximations are inevitable. Some approximations are shown that ease the numerical burden significantly but still faithfully reproduce the scattering rates expected from Fermi's golden rule. It is worth mentioning that all presented approximations have been carefully assessed. Approximations that violate Pauli blocking, underestimate effects of nonlocality of quantum mechanics or miss important characteristics of the scattering potentials have been avoided. Coherent effects or the balance between incoherent and coherent QCL physics are to the best of the authors' knowledge not affected by the presented approximations. A detailed discussion of the validity of approximations in NEGF can be found in Ref. 120. This section ends with a brief discussion on general approximations of scattering self-energies and some numerical details for efficient implementations.

#### 1. Scattering from longitudinal acoustic phonons

To avoid many-particle Green's functions, the phonon gas is assumed to remain unchanged by the electron propagation. This requires applying the perturbation potential of Eq. (66) in second order and to group the phonon creation with annihilation operators.^{216} The resulting products of phonon and electron operators are then translated into the electron and phonon Green's functions in Eq. (150) for the lesser and retarded self-energies^{195,209,211}