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.

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/AlSb4,5 and aluminum-free InGaAs/GaAsSb QCLs operating in the mid-infrared6 and terahertz7 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.

FIG. 1.

Conduction band profile and probability densities for a mid-infrared QCL9 lasing at 5 μm. Only the relevant energy states are displayed. The upper and lower laser levels are indicated by bold solid lines. Furthermore, the depopulation (dashed line) and injector (thin solid line) state probability densities are shown. The rectangle denotes a single QCL period.

FIG. 1.

Conduction band profile and probability densities for a mid-infrared QCL9 lasing at 5 μm. Only the relevant energy states are displayed. The upper and lower laser levels are indicated by bold solid lines. Furthermore, the depopulation (dashed line) and injector (thin solid line) state probability densities are shown. The rectangle denotes a single QCL period.

Close modal

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 theory34 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.

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.

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

0={22[1m(z)(x2+y2)+z1m(z)z]+V(z)E}ψ3D(x,y,z),
(1)

where ψ3D and E denote the wave function and eigenenergy, respectively. Furthermore, m 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., |ψ3D|2dxdydz=1. Since V and the effective masses only depend on the z coordinate, we can make the ansatz

ψ3D(x,y,z)=S1/2ψk(z)exp(ikxx+ikyy).
(2)

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 |ψk|2dz=1. Insertion of Eq. (2) in Eq. (1) yields the Ben Daniel-Duke model

{22kx2+ky2m(z)22z1m(z)z+V(z)Ek}ψk(z)=0,
(3)

where the wave function ψk(z) and energy Ek 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, yielding the one-dimensional (1D) Schrödinger equation in its usual form

[22z1m(z)z+V(z)E]ψ(z)=0.
(4)

The eigenenergy E due to 1D electron confinement in z direction is related to the total energy Ek by Ek = E + Ekin, where

Ekin=2(kx2+ky2)/(2m)=2k2/(2m)
(5)

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 composition37 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 mm.39 

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=E0iγq/2. The meaning of the imaginary term iγq/2 for quasi-stationary states can be understood by considering that the time dependence of a stationary wave function ψ is given by exp(iEt/). Consequently, for a complex eigenenergy E the density probability |ψ|2 of the electron in the quantum well system decays as exp(γqt), thus γq is the probability density decay rate of the quasi-bound state, and τq=γq1 can be associated with the lifetime of the particle in the quantum well.40 

FIG. 2.

Quasi-bound states in the conduction band of a QCL heterostructure.

FIG. 2.

Quasi-bound states in the conduction band of a QCL heterostructure.

Close modal

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 τ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 ψ=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 Ep, 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 Lp, and corresponding shifts in energy. For example, the solutions in the right-neighboring period are given by ψi(z)=ψi(zLp),Ei=EiEp, where ψi and Ei are the wave function and eigenenergy of the ith solution of Eq. (4) obtained in the central period.

FIG. 3.

Finite simulation window with artificial boundary conditions ψ=0. The central QCL period of the simulated structure is indicated, as well as the energy interval used for determining the bound states, corresponding to an energy period. The states in the other periods are found by adequate shifts in position and energy.

FIG. 3.

Finite simulation window with artificial boundary conditions ψ=0. The central QCL period of the simulated structure is indicated, as well as the energy interval used for determining the bound states, corresponding to an energy period. The states in the other periods are found by adequate shifts in position and energy.

Close modal

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

Ẽi=EiV|ψi|2dz
(6)

is a useful quantity.43 An automated selection of the relevant subbands can then be implemented by choosing the states with the lowest energies Ẽi in each period, considering all subbands which contribute significantly to the carrier transport.

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

m(E)=m2αE[1(14αE)1/2],
(7a)
m(E)=m[1+(2α+β)E].
(7b)

Here, α and β are nonparabolicity parameters, which are defined in the framework of a 14-band k·p calculation.51,52 Approximately, we have α=(Eg+Δso/3)1,49 where Eg and Δ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(z)=EV(z). Since the semiconductor material changes along the growth direction, also m,α, and β are z dependent. This approach is valid for moderate energies E1/α. The fact that Eq. (7a) even breaks down for αE>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 can assume relatively large values (see Fig. 3). This issue can be avoided by using Eq. (7a) for E<0, and its second order Taylor expansion

m(E)=m(1+αE),
(8)

for E>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 

FIG. 4.

Band structure of gallium arsenide (GaAs), obtained with NEMO5 in 20 band empirical tight binding.44 Shown are the valence bands (dashed lines), the conduction band (solid line), and the parabolic dispersion relation used for the Γ valley (dotted line).

FIG. 4.

Band structure of gallium arsenide (GaAs), obtained with NEMO5 in 20 band empirical tight binding.44 Shown are the valence bands (dashed lines), the conduction band (solid line), and the parabolic dispersion relation used for the Γ valley (dotted line).

Close modal

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(E), which is however z dependent. This problem can be avoided by defining an average effective mass for the ith subband

mi=m(z){1+[2α(z)+β(z)][EiV(z)]}|ψi(z)|2dz,
(9)

where the wave function is assumed to be normalized, |ψi|2dz=1, and Ei denotes the corresponding subband eigenenergy. This lowest order implementation of nonparabolicity thus yields different, albeit constant in-plane masses for each subband.

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 method41,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.

FIG. 5.

Numerical solution of the one-dimensional effective mass Schrödinger equation: (a) Transfer matrix approach and (b) finite difference method.

FIG. 5.

Numerical solution of the one-dimensional effective mass Schrödinger equation: (a) Transfer matrix approach and (b) finite difference method.

Close modal

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 znz<zn+Δn=zn+1 are approximated by constant values, e.g., Vn=V(zn),mn=m(zn), resulting in a jump VnVn+1,mnmn+1 at the border between the segments n and n + 1.41 Nonparabolicity can straightforwardly be implemented by using Eq. (7a) for m*; then mn=mn(E) depends on the eigenenergy E. The solution of Eq. (4) in segment n is given by

ψ(z)=Anexp[ikn(zzn)]+Bnexp[ikn(zzn)].
(10)

Here, kn=2mn(EVn)/ is the wave number (for E<Vn, we obtain kn=iκn=i2mn(VnE)/).41 The matching conditions for the wave function at a potential or effective mass discontinuity read

ψ(z0+)=ψ(z0),[zψ(z0+)]/m(z0+)=[zψ(z0)]/m(z0),
(11)

where z0+ and z0– denote the positions directly to the right and left of the discontinuity.61 Using Eqs. (10) and (11), the amplitudes An+1 and Bn+1 can be related to An and Bn by

(An+1Bn+1)=Tn,n+1(AnBn),
(12)

where the transfer matrix is with k̃n=kn/mn given by

Tn,n+1=Tnn+1Tn(Δn)=(k̃n+1+k̃n2k̃n+1eiknΔnk̃n+1k̃n2k̃n+1eiknΔnk̃n+1k̃n2k̃n+1eiknΔnk̃n+1+k̃n2k̃n+1eiknΔn).
(13)

Here,

Tn(Δn)=(eiknΔn00eiknΔn)
(14)

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

Tnn+1=12k̃n+1(k̃n+1+k̃nk̃n+1k̃nk̃n+1k̃nk̃n+1+k̃n)
(15)

is the potential step matrix for the interface between the segments n and n + 1 at z0 = zn+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, A0, B0 and AN, BN, can be obtained by multiplying the matrices for all segments

(ANBN)=TN1,NTN2,N1T0,1(A0B0)=(T11T12T21T22)(A0B0).
(16)

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 ψ=0 at the boundaries of our simulation window, corresponding to A0 + B0 = 0 and AN + BN = 0. The left boundary condition can, for example, be enforced by setting A0 = 1, B0 = −1. The right boundary condition AN + BN = 0 can then only be satisfied if the energy dependent matrix elements T11(E), T12(E), T21(E), and T22(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, ψN(E)=AN(E)+BN(E), and the eigenenergies are given by ψ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ΔE with sufficiently close spacing ΔE. The eigenenergies are located in intervals Em..Em+1 with ψN(Em)ψ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

Tn,n+1=Tn+1(Δn2)Tnn+1Tn(Δn2)=(k̃n+1+k̃n2k̃n+1eikn+Δnk̃n+1k̃n2k̃n+1eiknΔnk̃n+1k̃n2k̃n+1eiknΔnk̃n+1+k̃n2k̃n+1eikn+Δn),
(17)

where kn±=(kn±kn+1)/2 and kn=2mn(EVn)/,k̃n=kn/mn.43 Here, the band edge discontinuities occurring between wells and barriers must be treated separately by using the corresponding transfer matrix Eq. (15).

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 Δz is introduced, and the wave function ψ(z), potential V(z), and effective mass m(z) are represented by the corresponding values ψn, Vn, and mn on the grid points zn. First order derivatives are approximated by zψ(zn+Δz/2)Δψn+1/2=(ψn+1ψn)/Δz. Consequently, the term z(m)1zψ in Eq. (4) at z = zn can be expressed as (Δψn+1/2/mn+1/2Δψn1/2/mn1/2)/Δz. Using linear interpolation mn+1/2=(mn+mn+1)/2, we obtain the discretized form of Eq. (4)58,59

snψn1+dnψnsn+1ψn+1=Eψn,
(18)

with

sn=2Δz2(mn1+mn),
(19a)
dn=2Δz2(1mn1+mn+1mn+mn+1)+Vn.
(19b)
Here, the grid should be chosen so that the band edge discontinuities occur halfway between two adjacent grid points,58 as illustrated in Fig. 5(b). We again assume that the wave function is zero at the boundaries of our simulation window, ψ0=ψN=0. Writing Eq. (18) in matrix form yields an eigenvalue equation of the form (HEI)ψ=0, where ψ is the wave function vector with ψ=[ψ1,ψ2,,ψN1]T, I represents the identity matrix of size N – 1, and H is the Hamiltonian matrix with the non-zero elements Hn,n = dn, Hn,n–1 = −sn, Hn,n+1 = −sn+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 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.

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

e1z[ϵ(z)zṼ(z)]=e[nD(z)inis|ψi(z)|2].
(20)

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), ϵ(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 ψi(z). This charge distribution in the structure gives rise to space charge effects, resulting in an additional electrostatic potential energy Ṽ(z) which causes conduction band bending.67 

The total potential V in Eq. (4) is then given by V=V0+Ṽ. Here, V0=VcEpz/Lp, where Vc is the unbiased conduction band profile due to the varying material composition, thus describing the wells and barriers, and the term –Epz/Lp results from the applied bias. Since the energy drop across a period is given by the external bias, V(z0)V(z0+Lp)=Ep, we have Ṽ(z0)=Ṽ(z0+Lp). Due to the charge neutrality in each period, z0z0+Lpρdz=0, we furthermore obtain zṼ(z0)=zṼ(z0+Lp), i.e., Ṽ has the periodicity of Vc. Thus, we can restrict the solution of Eq. (20) to a single QCL period zϵ[z0,z0+Lp] and assume the boundary conditions Ṽ(z0)=Ṽ(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

s̃nṼn1d̃nṼn+s̃n+1Ṽn+1=ρn,
(21)

with

ρn=e[nD,ninis|ψi,n|2]
(22)

and

s̃n=12eΔz2(ϵn1+ϵn),d̃n=12eΔz2(ϵn1+2ϵn+ϵn+1).
(23)

Equation (21) is then solved over a single QCL period zϵ[z0,z0+Lp], with the grid points zn, n=0P which should coincide with the grid used for solving the Schrödinger equation Eq. (4). Applying the boundary conditions Ṽ0=ṼP=0, Eq. (21) can be written as a matrix equation MṼ=ρ, where Ṽ and ρ represent vectors with elements Ṽn and ρn, respectively, with n=1(P1). M is a matrix with the non-zero elements Mn,n=d̃n,Mn,n1=s̃n,Mn,n+1=s̃n+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

nis=miπ2kBTLln{1+exp[(μẼi)/(kBTL)]}.
(24)

Here, we assume that the electron distribution is described by the lattice temperature TL. Furthermore, μ is the chemical potential, kB denotes the Boltzmann constant, and mi is the effective mass associated with the ith subband, see Eq. (9). If nonparabolicity effects are neglected, mi can often be approximated by the value of the well material. In Eq. (24), the energy Ẽi defined in Eq. (6) is used instead of the eigenenergy Ei, 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.,

ns=z0z0+LpnDdz=inis.
(25)

This is done recursively by first determining a lower and an upper boundary value for μ where inis<ns and inis>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+Ṽ, where Ṽ has to be obtained by solving the Poisson equation (20). On the other hand, the wave functions ψ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 Ṽ=0, until the results for Ṽ,ψi, and Ei 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 QCL69 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 Ṽ will be similar as for the simplified model based on Eq. (24).

FIG. 6.

Conduction band profile and probability densities for a bound-to-continuum terahertz QCL.69 Shown are the results without considering space charge effects (solid lines) and for space charge effects taken into account, assuming thermally occupied subbands (dashed lines).

FIG. 6.

Conduction band profile and probability densities for a bound-to-continuum terahertz QCL.69 Shown are the results without considering space charge effects (solid lines) and for space charge effects taken into account, assuming thermally occupied subbands (dashed lines).

Close modal

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 

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

FIG. 7.

Typical waveguide resonator geometry of a QCL.

FIG. 7.

Typical waveguide resonator geometry of a QCL.

Close modal

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 exp(iωt), where ω denotes the angular frequency. The fields are then expressed as H(x,t)={Ĥ(x)exp(iωt)}, E(x,t)={Ê(x)exp(iωt)}, where Ĥ and Ê denote the complex amplitude vectors. All equations can easily be converted to the engineering convention [expressing the time dependence as exp(jωt)] by making the formal substitution ij. 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 ϵr(x). Under these conditions, Maxwell's equations simplify to

×Ĥ=iωϵ0ϵrÊ,
(26a)
×Ê=iωμ0Ĥ.
(26b)
B=μ0(H)=0 is then automatically fulfilled, as can be seen by taking the divergence of Eq. (26b). The equation D=(ϵ0ϵrE)=ρ delivers the charge density ρ and is not needed to compute 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 element73 and finite difference time domain35 method, which also have been applied to the simulation of QCL cavities.71,74

By eliminating Ê from Eqs. (26), a wave equation for Ĥ can be derived.75 In many cases, a waveguide geometry is used which does not depend on the longitudinal x direction, i.e., ϵr=ϵr(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)={Ĥy,z(y,z)exp(iβxiω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

(p2+q2)Ĥp+qϵrϵr(pĤqqĤp)=(β2k02ϵr)Ĥp,
(27)

with k0=ωμ0ϵ0=ω/c,p=y,z, and q = z,y.75 The longitudinal component Hx is obtained from H=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, constitutes an eigenvalue problem, and the corresponding solutions Ĥy,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., Hx = 0.

A waveguide mode is frequently characterized using three parameters, the overlap (or field confinement) factor Γ, the waveguide loss coefficient aw, and the mirror or outcoupling loss coefficient am, with the total power loss coefficient a = aw + am. Based on these quantities, the threshold gain is given by gth=a/Γ. 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 Γ, aw, and am 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{β}. 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 Ez component, the overlap factor is defined as78 

Γ=Sg|Êz|2dydz|Ê|2dydz,
(28)

where we integrate over the gain medium cross section area Sg 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

R=|neff1|2/|neff+1|2,
(29)

with the effective refractive index defined as neff=β/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 am is obtained from R=exp(amL) with the resonator length L, yielding

am=ln(R)/L.
(30)

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

If the waveguide width in lateral y direction significantly exceeds its thickness, the waveguide calculations can be reduced to a 1D problem with ϵr=ϵr(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 (25200μm).71 For TM modes, in addition to Hx = 0 we have approximately Hz0, and the y component is given by Hy(x,z)={Ĥy(z)exp(iβxiωt)}, where β denotes the complex propagation constant. The one-dimensional wave equation is then obtained from Eqs. (26) or Eq. (27) as

ϵrz(ϵr1zĤy)=(β2ϵrk02)Ĥy.
(31)

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

Êz=βωϵ0ϵrĤy,
(32a)
Êx=iωϵ0ϵrzĤy.
(32b)

From Eq. (31) we see that both Hy and (1/ϵr)zHy must be continuous, corresponding to the continuity of the field components Hy and, with Eq. (32b), Ex parallel to the layers.

Equation (31) can be brought into the same form as the one-dimensional Schrödinger equation (4) by associating Ĥy(z),ϵr(z), and ϵr1β2k02 with ψ(z),2m(z)/2, and EV(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, znz<zn+Δn=zn+1, with constant relative permittivities ϵr(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 ϵr,eff1=(Δbϵr,b1+Δwϵr,w1)/(Δb+Δw), where Δb and Δw denote the total thickness of the barriers and wells in the gain medium, respectively, and ϵr,b,ϵr,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 as83 

Ĥy(n)=Anexp[ikn(zzn)]+Bnexp[ikn(zzn)],
(33)

with kn=(ϵr(n)k02β2)1/2. The propagation through the segment n is then described by the matrix Tn(Δn) defined in Eq. (14). As mentioned above, Eq. (31) implies the continuity of Hy and (1/ϵr)zHy across layer boundaries. These matching conditions between two layers can be expressed in terms of matrix Eq. (15) by choosing k̃n=kn/ϵr(n). There is, however, one fundamental difference between Eqs. (4) and (31): Equation (31) generally has complex eigenvalues β2 since ϵr 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 Ĥy(z). For TM polarization, we obtain from the boundary value method79 

T=[{β}{ϵr(z)}+{β}{ϵr(z)}|ϵr(z)|2|Ĥy(z)|2dz]1×8π|β|2k0k0k01kz2/k02|Φ(kz)|2|Φ(kz)|2|k01kz2/k02Φ(kz)+βΦ(kz)|2dkz,
(34)

with the Fourier transforms

Φ(kz)=12πĤy(z)exp(ikzz)dz,Φ(kz)=12πĤy(z)ϵr(z)exp(ikzz)dz.
(35)

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 

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 

FIG. 8.

Illustration of various carrier transport models with different levels of complexity and accuracy: (a) Rate equation approach; (b) ensemble Monte Carlo method; (c) 3D density matrix/NEGF approach.

FIG. 8.

Illustration of various carrier transport models with different levels of complexity and accuracy: (a) Rate equation approach; (b) ensemble Monte Carlo method; (c) 3D density matrix/NEGF approach.

Close modal

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 form114,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.

Table I.

Classification of carrier transport modeling techniques. The corresponding section number is given in brackets.

  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.

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 

The rate equations for a laser are given by88 

dtnis=jiτji1njsτi1nis+ji(Wijoptnis+Wjioptnjs).
(36)

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 τji1, where τi1=jiτij1 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 |ωij|=|EiEj|/ 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 formation127 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

dtnis=jiτ̂ji1njsτi1nis+ji(Ŵijoptnis+Ŵjioptnjs),
(37)

with i,j=1..N, where N is the number of subbands in each period. Here, τ̂ji1=nτj,i+nN1 includes the transitions to all equivalent levels in the different periods, and analogously for Ŵjiopt. The total sheet density in each period is determined by the doping sheet density

ns=i=1Nnis.
(38)

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, Ŵijopt=Ŵjiopt=0. The subband populations nis(i=1..N) can then be found by solving the linear equation system Eq. (37) with dtnis=0 and i=1..(N1), 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 Ŵijopt 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 η2 of the current density J is injected into the upper laser level, and η1 into the lower laser level. Furthermore, only the transition rates τ211,τ201, and τ101 are included, with the inverse upper and lower laser level lifetimes τ21=τ211+τ201,τ11=τ101. With Eq. (36), the rate equations for the upper and lower laser level are then obtained as

dtn2s=η2J/eτ21n2sσ21optI(n2sn1s),dtn1s=η1J/eτ11n1s+τ211n2s+σ21optI(n2sn1s),
(39)

where σ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 ω21=(E2E1)/. If lasing is neglected, I = 0, the steady state population inversion is with dt = 0 in Eq. (39) obtained as89 

n2sn1s=Je[η2τ2(1τ1τ21)η1τ1].
(40)

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 τ01 and τ12, respectively.90 

FIG. 9.

Three-level system model for a QCL.

FIG. 9.

Three-level system model for a QCL.

Close modal

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 by95 

tρ21=iω21ρ21i1d12EzΔ21γ21ρ21,tΔ21=2i1(d12ρ21d12ρ21)EzγE(Δ21Δ21eq),
(41)

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=e1|z|2 is the corresponding dipole matrix element of the laser transition. Furthermore, ρij(x,t)=ρji(x,t)=i|ρ̂(x,t)|j are the density matrix elements, and Δ21(x,t)=ρ22(x,t)ρ11(x,t) is the inversion. The density matrix can be normalized so that ρii=nis/ns gives the relative population of subband i. Dissipative processes are phenomenologically included by adding decay terms with relaxation rates γE and γ21, describing the energy relaxation and dephasing, respectively. Δ21eq is the equilibrium inversion which the system approaches for Ez = 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 d12, γE,γ21,ω21, and Δ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

(x2c2n02t2)Ez=ϵ01c2t2Pz.
(42)

Pz(x,t) is the polarization component in z direction due to the lasing transition, given by Pz=(ns/Lp)(d12ρ21+d12ρ21). Here, Lp is the length of a single QCL period and ns /Lp thus corresponds to the average electron concentration.97 Furthermore, n0 denotes the refractive index of the gain medium material.

Ez and ρ21 are typically expressed by their slowly varying envelope functions

Ez(x,t)=12Êz(x,t)exp[i(k21xω21t)]+c.c.,ρ21(x,t)=η21(x,t)exp[i(k21xω21t)],
(43)

where k21=ω21n0/c and c.c. denotes the complex conjugate. Inserting Eq. (43) into Eq. (41) and neglecting the rapidly oscillating terms exp(±2iω21t), we obtain the density matrix equations in the rotating wave approximation95,

tη21=i(2)1d12ÊzΔ21γ21η21,
(44a)
tΔ21=i1(η21d12Êzη21d12Êz)γE(Δ21Δ21eq).
(44b)
Inserting Eq. (43) into Eq. (42) yields with the slowly varying amplitude approximation |x2Êz||k21xÊz|,|t2Êz||ω21tÊz|95,97

xÊz+n0c1tÊz=iω21Γϵ0cn0nsLpd12η2112aÊz,
(45)

furthermore assuming |t2η21|,|ω21tη21||ω212η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, τ10, resulting in n1s=0 in Eq. (39).97 The inversion is then directly given by the upper laser level population, Δ21=n2s. By comparison of Eq. (44) with Eq. (39), we obtain Δ21eq=τ2η2J/(nse),γE=τ21.

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 tη21=0

η21=i2γ21d12ÊzΔ21.
(46)

Multiplying Eq. (45) from left with Êz and adding the complex conjugate, we obtain with Eq. (46)

xI+n0c1tI=ΓgIaI,
(47)

with the power gain coefficient

g=ω21γ21ϵ0cn0nsLp|d12|2(ρ22ρ11).
(48)

Here, we have replaced the electric field by the optical intensity I=ϵ0cn0|Êz|2/2. The optical power inside the resonator is given by P=ISg/Γ, where Sg is the cross section area of the gain medium and Sg/Γ 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 am, Eq. (30). Equation (47) then simplifies to

n0c1tI=ΓgIaI.
(49)

The transition rate due to the optical field is with tΔ21|opt=2tρ22|opt=2tρ11|opt obtained from Eqs. (44b) and (46) as

tρ22|opt=12i1(η21d12Êzη21d12Êz)=12γ21ϵ0cn0|d12|2(ρ22ρ11)I.
(50)

The contribution ρ22I leading to a reduction of ρ22 is due to stimulated emission, while the contribution ρ11I corresponds to absorption. For a slightly detuned optical field at a frequency ωω21, Eqs. (48) and (50) can be adapted by replacing γ211 with πL(ω), where L(ω) is the Lorentzian lineshape function95 

L(ω)=1πγ21γ212+(ω|ω21|)2.
(51)

Thus, we obtain with the sheet densities n1s=nsρ11, n2s=nsρ22,

g=πωϵ0cn0Lp|d12|2(n2sn1s)L(ω),
(52)

and

tn2s|opt=πϵ0cn02|d12|2(n2sn1s)IL(ω).
(53)

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

Wijopt=Wjiopt=πϵ0cn02|d12|2IL(ω).
(54)

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 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, Vcos(ω0t), the carrier energy is changed by ω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 

FIG. 10.

Different classes of scattering processes: (a) elastic scattering; (b) inelastic scattering; (c) carrier-carrier scattering.

FIG. 10.

Different classes of scattering processes: (a) elastic scattering; (b) inelastic scattering; (c) carrier-carrier scattering.

Close modal
FIG. 11.

Simulation results for the spectral gain vs. frequency, as obtained by evaluating electron-phonon (e-p), electron-impurity (e-i), and electron-electron (e-e) scattering (solid curves) and neglecting e-e (dashed curves) or e-i (dotted curves) scattering. (a) 3.4 THz resonant phonon depopulation structure; (b) 3.5 THz bound-to-continuum structure. Reprinted with permission from J. Appl. Phys. 105, 123102 (2009). Copyright 2009 AIP Publishing LLC.107 

FIG. 11.

Simulation results for the spectral gain vs. frequency, as obtained by evaluating electron-phonon (e-p), electron-impurity (e-i), and electron-electron (e-e) scattering (solid curves) and neglecting e-e (dashed curves) or e-i (dotted curves) scattering. (a) 3.4 THz resonant phonon depopulation structure; (b) 3.5 THz bound-to-continuum structure. Reprinted with permission from J. Appl. Phys. 105, 123102 (2009). Copyright 2009 AIP Publishing LLC.107 

Close modal

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 to a final state |jk in the QCL heterostructure, where k=(kx,ky)T and k=(kx,ky)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 ψ3D,i=S1/2ψi(z)exp(ikr) with r=(x,y)T and Eik=Ei+2k2/(2mi) where ψi(z) and Ei 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

Vjk,ik=jk|V|ik=S1SVψjψiexp[i(kk)r]d2rdz.
(55)

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

(Vjk,ik+(Q)Vjk,ik(Q))=jk|(V0exp(iQx)V0exp(iQx))|ik=S1Sψjψi(V0exp(iQx)V0exp(iQx))×exp[i(kk)r]d2rdz.
(56)

Assuming bound wave functions ψi,j(z) in Eqs. (55) and (56), the integration over the z coordinate can be taken from to . 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 to also in x and y direction.

The transition rate from an initial state |ik to a final state |jk is obtained from Fermi's golden rule,61 given by

Wik,jk=2π|Vjk,ik|2δ(EjkEik)
(57)

for elastic scattering processes and

Wik,jk±(Q)=2π|Vjk,ik(Q)|2δ(EjkEik±ω0)
(58)

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

|k|=k=(mjmik2+2mjEiEj2)1/2.
(59)

Analogously, energy conservation yields for inelastic scattering

k=(mjmik2+2mjEiEjω02)1/2.
(60)

The computation of the total transition rates from an initial state |ik to a final subband j involves the summation over wave vectors. These sums can be converted to integrals, introducing a factor of Ld/(2π) per dimension where the device length Ld 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 which is two-dimensional. It can furthermore be advantageous to express k in polar coordinates |k| and ϕ, and introduce a kinetic energy variable Ejkin=2|k|2/(2mj). Thus, we obtain

kS(2π)2d2k=Smj(2π)2002πdϕdEjkin.
(61)

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 Ωc be approximated by

QΩc(2π)3d3Q.
(62)

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ωQt), with the amplitude U, wave vector Q, and angular frequency ωQ. The normal modes are the solutions of u for which the lattice uniformly oscillates at a single frequency ωQ, and a phonon corresponds to an elementary vibrational motion. The associated relation between wave vector Q and frequency ω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 ωQ=0 for Q = 0. For optical modes, two consecutive atoms in the same unit cell move in opposite direction, and ωQ at Q = 0 typically corresponds to infrared optical frequencies. Furthermore, the wave propagation can be predominantly longitudinal (QU) or transverse (QU); 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 Nu – 1 longitudinal optical (LO) and 2(Nu1) transverse optical (TO) branches, where Nu denotes the number of atoms per unit cell. For example, GaAs has two atoms per unit cell, i.e., Nu = 2.

FIG. 12.

Schematic illustration of the dispersion relation for an acoustic and optical phonon branch. This simplified example corresponds to a one-dimensional lattice with Nu = 2 and lattice constant a.

FIG. 12.

Schematic illustration of the dispersion relation for an acoustic and optical phonon branch. This simplified example corresponds to a one-dimensional lattice with Nu = 2 and lattice constant a.

Close modal

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 Nu2) 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 Nu2) 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ωQt)

E=Nam|tu|2t=12NamU2ωQ2,
(63)

where the number of atoms with mass m in a crystal lattice of volume Ωc and density ρc is Na=Ωcρc/m. Thus, we obtain the amplitude for a mode occupied by a single phonon of energy E=ωQ

U=2ΩcρcωQ.
(64)
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 Shockley138 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=δVc. The corresponding constant of proportionality is the deformation potential Ξ (for valleys without rotational symmetry, Ξ becomes a tensor), i.e., V=Ξu for small displacements.61 We obtain

V=ΞQUcos(QxωQt),
(65)

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

V=M(Q){exp[i(QxωQt)]+exp[i(QxωQt)]},
(66)

with the amplitude

M(Q)=ΞQ22ΩcρcωQ=ΞQ2Ωcρcvs.
(67)

Equation (58) gives the transition rate for an electron in an initial state |ik, which is scattered to a final state j with a wave vector k 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

Vjk,ik±=M(Q)SSexp[i(k±qk)r]d2rFji±(qz),
(68)

with the form factor

Fji±(qz)=ψj(z)exp(±iqzz)ψi(z)dz.
(69)

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

Wik,jk±(Q)=(2π)3S(NQ+12±12)|M(Q)|2|Fji(qz)|2×δ(kqk)δ(EjkEik±ωQ),
(70)

where Wik,jk+ and Wik,jk refer to emission and absorption, respectively. Here, we have used that |Sexp[i(k±qk)r]d2r|2 can be approximated by 4π2Sδ(kqk) 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

NQ=[exp(ωQkBTL)1]1.
(71)

For emission, the factor NQ + 1 is used to also include spontaneous emission processes.

The total transition rate from a given initial state |ik to a subband j is obtained by summing over all wave vectors k 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 to a subband j132 

Wik,j±=Ξ28π2ρcvsQ(NQ+12±12)|Fji(qz)|2×δ[Ej,kqEik±ωQ]d3Q.
(72)

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 ±ωQ in Eqs. (70) and (72). Besides, we can approximate Eq. (71) as

NQNQ+1kBTLωQkBTLvsQ,
(73)

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

Wik,jk±(Q)=Ξ24π3kBTLΩcρcvs2S|Fji(qz)|2×δ(kqk)δ(EjkEik).
(74)

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

Wik,j±=Ξ2kBTLmj4πρcvs23|Fji(qz)|2dqz.
(75)

The final wave vector magnitude k is then given by Eq. (59). With Eq. (69), we obtain

|Fji(qz)|2dqz=2π|ψj(z)|2|ψi(z)|2dz.
(76)

Thus, the total transition rate is given by

Wik,j=Ξ2kBTLmjρcvs23|ψj(z)|2|ψi(z)|2dz,
(77)

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, δVc=DTOu=DTOu. Here, DTO is the optical deformation potential field which assumes a different value for different kinds of intra- and intervalley transitions, and DTO is the component parallel to u. For optical phonons, ωQ reaches an extremum ωQ=ωTO for Q = 0, and for small Q, the dispersion relation can thus be approximated by ωQ=ωTO (see Fig. 12). The scattering rate for an electron from an initial state |ik to a final state |jk due to phonons with wave vector Q is again given by Eq. (70), with

M(Q)=DTO2ΩcρcωTO.
(78)

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

Wik,j±=DTO28π2ρcωTO(NQ+12±12)|Fji(qz)|2×δ[Ej,kqEik±ωTO]d3Q.
(79)

Approximating ωQ=ωTO, the phonon number NQ in Eq. (79) does not depend on Q, and thus we can set NQ = NPh. With Eq. (76), the result then simplifies to

Wik,j±=DTO2mj2ρcωTO2(NPh+12±12)×|ψj(z)|2|ψi(z)|2dz,
(80)

for EjEikω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 ωQ=ωLO for small Q (see Fig. 12), where ω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 mA and mB in a unit cell, these effective charges are given by ±qeff with qeff2=ϵ0ωLO2μnc1(ϵr,1ϵr,01).61 Here, μ1=mA1+mB1 is the inverse reduced mass, and nc is the number of cells per unit volume. Furthermore, ϵr,0 denotes the static dielectric constant and ϵr, 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ω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 ρc is replaced by μnc. The resulting ionic polarization is Pion=pnc, and the electric field E can be obtained from the material equation D=ϵ0E+Pion with D = 0 due to the absence of free carriers. With the potential V=eEdx, the resulting scattering rate from a state |ik to |jk is again given by Eq. (70), where

M(Q)=12Qeϵ01(ϵr,1ϵr,01)2ωLOΩc.
(81)

The total transition rate from a given initial state |ik to a subband j is again obtained by summing over all wave vectors k and Q, using Eqs. (61) and (62). For EjEikωLO, we obtain

Wik,j±=mjωLOe28π22ϵ0(ϵr,1ϵr,01)(NPh+12±12)02πJ(q)dθ,
(82)

and Wik,j±=0 otherwise. Here,

J(q)=πqψj(z)ψj(z)ψi(z)ψi(z)×exp(q|zz|)dzdz
(83)

with the in-plane phonon wave vector magnitude

q=(k2+k22kkcosθ)1/2,
(84)

where k is given by Eq. (60) with ω0=ω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 ϵr,0, the use of bulk modes is a good approximation.101 

Screening can be included by changing Eq. (81) to141,142

M(Q)=Q2(Q2+qs2)eϵ01(ϵr,1ϵr,01)2ωLOΩc,
(85)

where qs is the inverse screening length. For Eq. (83), we then obtain with q̃=q2+qs2

J(q)=πq̃ψj(z)ψj(z)ψi(z)ψi(z)exp(q̃|zz|)×(1|zz|qs22q̃qs22q̃2)dzdz.
(86)

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

qs=nee2ϵ0ϵr,0kBTL,
(87)

where ne is the averaged three-dimensional electron density. However, Debye screening is generally only valid for high temperatures,143 since qs diverges for TL0.

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.

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 scatters to a final state |jk, accompanied by a transition of a second electron from a state |i0k0 to |j0k0. 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 to a subband j is then obtained as141,150

Wik,j=m4π3Si0,j0,k0fi0(k0)02πdθ|Mii0jj0(q)|2,
(88)

with the in-plane cross section area S and carrier distribution function fi0(k0) for the state |i0k0. θ is the angle between g=k0k and g=k0k, and q=kk (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 ϵ=ϵ0ϵr,0 given by

Vii0jj0b(q)=e22ϵqdzdz[ψi(z)ψi0(z)×ψj(z)ψj0(z)exp(q|zz|)].
(89)

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 qs. For the transition matrix element Mii0jj0, we then obtain

|Mii0jj0(q)|2=12|qq+qsVii0jj0b(q)|2,
(90)

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(Ei+Ei0EjEj0)/2

q=12[2g2+g022g(g2+g02)1/2cosθ]1/2.
(91)

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 then110,152

qs=e22ϵmπ2fi(k=0).
(92)

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 system154 

Vii0jj0s=Vii0jj0b+mnVimjnbΠmnVmi0nj0s.
(93)

Here, Πmn(q) is the polarizability tensor, given in the long wavelength limit (q0) by

Πmn={nmsnnsEmEn,mn,mπ2fn(0),m=n.
(94)

Commonly, simplified screening models using a constant screening wave number qs are employed to avoid the numerical load associated with solving Eq. (93).110,152Vii0jj0s is then obtained by replacing the prefactor e2/(2ϵq) in Eq. (89) by e2/[2ϵ(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 model152 with

qs=e22ϵmπ2ifi(k=0),
(95)

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 Vii0jj0sand the “exchange” matrix element Vii0j0js.153 Accounting for this exchange effect, we obtain150,153

|Mii0jj0|2=pa2[|Vii0jj0s(q)|2+|Vii0j0js(q+)|2]+pp2|Vii0jj0s(q)Vii0j0js(q+)|2,
(96)

where

q±=12[2g2+g02±2g(g2+g02)1/2cosθ]1/2,
(97)

and pa = pp = 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 pa = 1/2, pp = 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 pa = 1, pp = 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 pa = pp = 1/2 (solid curve), pa = 1, pp = 0 (dotted curve), pa = 1/2, pp = 0 (dashed curve), and pa = pp = 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.

FIG. 13.

EMC simulation results for the spectral gain vs. frequency of a bound-to-continuum terahertz QCL69 at a lattice temperature TL = 10 K. (a) Comparison of different screening models: RPA (solid curve), modified single subband model for all matrix elements (dashed curve) or for the intrasubband elements only (dotted curve), thus treating intersubband elements as unscreened. (b) Comparison of different spin implementations. Shown are results for fully taking into account (solid curve) and ignoring the exchange effect (dotted curve), and ignoring parallel spin collisions (dashed curve). For comparison, the result obtained with no electron-electron scattering is also displayed (dashed-dotted curve). Reprinted with permission from J. Appl. Phys. 107, 013104 (2010). Copyright 2010 AIP Publishing LLC.68 

FIG. 13.

EMC simulation results for the spectral gain vs. frequency of a bound-to-continuum terahertz QCL69 at a lattice temperature TL = 10 K. (a) Comparison of different screening models: RPA (solid curve), modified single subband model for all matrix elements (dashed curve) or for the intrasubband elements only (dotted curve), thus treating intersubband elements as unscreened. (b) Comparison of different spin implementations. Shown are results for fully taking into account (solid curve) and ignoring the exchange effect (dotted curve), and ignoring parallel spin collisions (dashed curve). For comparison, the result obtained with no electron-electron scattering is also displayed (dashed-dotted curve). Reprinted with permission from J. Appl. Phys. 107, 013104 (2010). Copyright 2010 AIP Publishing LLC.68 

Close modal

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 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,y,z)T is with r=(x,y)T and ϵ=ϵ0ϵr,0 given by the Coulomb potential

V=(4πϵ)1e2(|rr|2+|zz|2)1/2.
(98)

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 as152 

Vijb(z,q)=e22ϵqψi(z)ψj(z)exp(q|zz|)dz.
(99)

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

Wik,jk=2πS|Vijb(z,q)|2nD(z)dzδ(EjkEik)=πe42Sϵ2q2Fij(q)δ(EjkEik),
(100)

with the form factor

Fij(q)=dznD(z)×[ψi(z)ψj(z)exp(q|zz|)dz]2.
(101)

Here, we have summed over the donors to include the effect of all ionized impurities, corresponding to an integration SnD(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, the final wave vector magnitude k is given by Eq. (59).

The total transition rate from a given initial state |ik to a subband j is found by summing over all final wave vectors k using Eq. (61). We obtain54,152

Wik,j=mje44πϵ230πFij(q)q2dθ,
(102)

with q(θ)=|kk|=(k2+k22kkcosθ)1/2. With Eq. (59), this yields

q(θ)=[(1+mjmi)k2+q022k(mjmik2+q02)1/2cosθ]1/2,
(103)

where q02=2mj(EiEj)/2. Values of θ resulting in complex q(θ) 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,q) for impurity scattering from the bare elements, Vijb(z,q) in Eq. (99). We obtain151,152

Vijs=Vijb+mnVimjnbΠmnVmns,
(104)

with Vimjnb and Πmn defined in Eqs. (89) and (94), respectively. In Eq. (100), Vijb(z,q) has then to be replaced by Vijs(z,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 qs 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 Herring158 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 Δ(x,y)=Δ(r) from its average position z0 as illustrated in Fig. 14. Since it is not feasible to measure Δ(r) directly, the interface roughness is characterized by its standard deviation Δ and correlation length Λ. Typically, a Gaussian autocorrelation function is assumed162 

Δ(r)Δ(r)=1SΔ(r)Δ(r+d)d2r=Δ2exp(d2Λ2),
(105)

with d=rr. A change Δ(r) of the interface position corresponds to a perturbing potential

V=±Vo[H(zz0)H(zz0Δ(r))],
(106)

where z0 and Vo 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 Δ(r) is small, i.e., on the order of a monolayer, we can approximate ψi,j(z0+Δ(r))ψi,j(z0). Thus, we obtain with Eq. (55)

Vjk,ik=±VoSψi(z0)ψj(z0)d2rΔ(r)exp(iqr),
(107)

and

|Vjk,ik|2=Vo2S2|ψi(z0)ψj(z0)|2×d2d[exp(iqd)d2rΔ(r)Δ(r+d)],
(108)

with q=kk. Assuming a Gaussian autocorrelation Eq. (105), we obtain from Eq. (57) the scattering rate

Wik,jk=2π2SVo2Δ2Λ2n|ψi(zn)ψj(zn)|2×exp(14Λ2q2)δ(EjkEik),
(109)

where we additionally sum over all interfaces, located at positions zn. Due to the energy conservation Eik=Ejk, the final wave vector magnitude k is given by Eq. (59).

FIG. 14.

Interface between two heterostructure layers consisting of different materials. The local deviation from a perfect interface due to roughness is expressed by Δ(x,y).

FIG. 14.

Interface between two heterostructure layers consisting of different materials. The local deviation from a perfect interface due to roughness is expressed by Δ(x,y).

Close modal

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

Wik,j=mj3Vo2Δ2Λ2n|ψi(zn)ψj(zn)|2×0πdθexp(14Λ2q2)H(q2),
(110)

where q(θ) 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 Λ10nm and Δ0.1nm.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 Δ=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.

FIG. 15.

Simulation results for the spectral gain vs. frequency, as obtained for a correlation length Λ=10nm and different values of the interface roughness mean height Δ. (a) 3.4 THz resonant phonon depopulation structure; (b) 3.5 THz bound-to-continuum structure. Reprinted with permission from J. Appl. Phys. 105, 123102 (2009). Copyright 2009 AIP Publishing LLC.107 

FIG. 15.

Simulation results for the spectral gain vs. frequency, as obtained for a correlation length Λ=10nm and different values of the interface roughness mean height Δ. (a) 3.4 THz resonant phonon depopulation structure; (b) 3.5 THz bound-to-continuum structure. Reprinted with permission from J. Appl. Phys. 105, 123102 (2009). Copyright 2009 AIP Publishing LLC.107 

Close modal

3. Alloy scattering

In a ternary semiconductor alloy AxB1–xC such as AlxGa1–xAs, 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 by36,168

Vc(x)=nVn(xxn),
(111)

where xn denotes the position of unit cell n. Furthermore, Vn = VAC (Vn = VBC) if the unit cell contains an atom of type A (B), where VAC and VBC 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,

Vc,VCA=n[xVAC(xxn)+(1x)VBC(xxn)].
(112)

The perturbing potential corresponds then to the difference between the real conduction band energy Vc(x) and Vc,VCA obtained by the virtual crystal approximation

V=VcVc,VCA=ncnδV(xxn),
(113)

with δV(x)=VAC(x)VBC(x) and cn=(1x)(cn=x) if the unit cell contains an atom of type A (B).

The scattering potential δ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 ψ(x) which are slowly varying on the atomic scale, δV can be modeled as a δ function potential36 

δV(x)=δVΩ0δ(x),
(114)

where Ω0 denotes the unit cell volume, given by Ω0=a3/4 for a zinc-blende structure with lattice constant a. Furthermore, δ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, δ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=kk

Vjk,ik=S1δVΩ0ncnψj(zn)ψi(zn)exp(iqrn),
(115)

and the square modulus is given by

|Vjk,ik|2=S2(δVΩ0)2nm{cncmexp[iq(rnrm)]×ψj(zn)ψj(zm)ψi(zn)ψi(zm)}.
(116)

For calculating the average of |Vjk,ik|2, the expectation value cncm 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 cncm=0 for mn and cncm=x(1x) for m = n. Furthermore replacing the sum over the unit cell positions xnf(zn) by an integral Ω01f(z)d3x=Ω01Sf(z)dz, we obtain

|Vjk,ik|2=S1Ω0(δV)2x(1x)×|ψj(z)ψi(z)|2dz.
(117)

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

Wik,jk=2πSδ(EjkEik)×Ω0(δV)2x(1x)|ψjψi|2dz,
(118)

where Ω0(z),δ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, the final wave vector magnitude k is given by Eq. (59). The total transition rate from a given initial state |ik to a subband j is with Eq. (61)

Wik,j=mj3Ω0(δV)2x(1x)|ψjψi|2dz
(119)

for Ej<Eik, and 0 otherwise.

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

fiFD(Eik)={exp[(EikEiF)/(kBTi)]+1}1,
(120)

where EiF is here a “quasi” Fermi energy describing the kinetic energy distribution of the electrons within the subband i,132 and Ti 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+2k2/(2mi) 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 obtain132 

(τij(m))1=EifiFD(E)Wik,j(m)dEEifiFD(E)dE.
(121)

Here, niE has been omitted in the enumerator and denominator because it is constant, niE=mi/(π2) for EEi. 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., τij1=m(τij(m))1.

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

nis=EifiFD(E)niEdE=mikBTiπ2{EiFEikBTi+ln[1+exp(EiFEikBTi)]}.
(122)

Furthermore expressing the carrier energy in terms of the in-plane wave vector, we obtain with dE=2(mi)1kdk132 

(τij(m))1=0fiFD[Ei+2k2/(2mi)]Wik,j(m)kdkπnis.
(123)

An additional factor {1fjFD[Ei+2k2/(2mi)±E0]} 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, E0 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)

EiFEi=kBTiln[exp(nisπ2mikBTi)1].
(124)

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

fiMB(Eik)=exp[(EikEiF)/(kBTi)].
(125)

Under this condition, Eq. (123) simplifies to

(τij(m))1=2mikBTi0exp(2k22mikBTi)Wik,j(m)kdk.
(126)

Numerically, the integration of Eq. (123) or Eq. (126) is performed from Ei up to a sufficiently large maximum value, e.g., the highest value of the simulated potential profile V(z).132 

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 solution92 

nis=jinjsτji1/jiτij1,
(127)

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, ξnis,new+(1ξ)nis,old, as input for the next iteration,91 where the relaxation parameter is typically chosen as ξ=0.5. The simulation has converged when the obtained sheet densities do not significantly change anymore between iterations.

3. Kinetic energy balance method

Since Ti 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, Ti = TL. However, the electron temperature can significantly exceed TL in quantum cascade lasers.174 Thus, the rate equation model has been extended to account for electron heating, typically assuming an identical electron temperature Ti = Te for all subbands.173 The kinetic electron energy generation rate per period and unit device in-plane cross section is given by173 

RE=i=1Njmnis(τij(m))1(EiEj+E0(m)),
(128)

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)=ELO) 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 Te where RE = 0 in Eq. (128). The rate equation is self-consistently solved as described above for an initial guess of Te, e.g., Te = TL. Based on the obtained nis and τij(m), RE in Eq. (128) is calculated. This procedure is repeated and the guess for Te is iteratively improved until RE0 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 Ti in the individual subbands.93 

Semiclassically, the carrier transport between states |ik in a quantum well system is given by the Boltzmann equation176 

dtfik=jk(Wjk,ikfjkWik,jkfik),
(129)

with the scattering rates Wik,jk=mWik,jk(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 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 ik and jk, 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 Ne104..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 tsim, which must be chosen long enough to ensure convergence to the stationary solution. The simulation is divided into subintervals Δt, where the scattering dynamics is subsequently evaluated. Δ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 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=2k2/(2mi) instead. The kinetic energy grid then divides the energy axis into segments n of widths ΔE(n) centered around discrete energies Ekin(n) with k(n)=(2miEkin(n))1/2/. 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̃ik,j(m) for time dependent scattering rates and compensating for the too high value by introducing artificial “self-scattering,” as described further below.

FIG. 16.

Schematic diagram of the EMC algorithm

FIG. 16.

Schematic diagram of the EMC algorithm

Close modal

An important quantity is the carrier distribution function, which can directly be obtained from Ni(n), denoting the number of simulated electrons in the nth 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/(π2),61 the number of available states in the nth energy cell is niEΔE(n)S. The simulated device in-plane cross section is with the sheet doping density per period ns given by S=Ne/(nsNp), where Np corresponds to the number of periods over which the Ne simulated electrons are distributed. The carrier distribution function is then approximately given as

fi(Ekin,t)=nsNpNi(n)(t)niEΔE(n)Ne,
(130)

where n indicates the energy cell containing the value Ekin. Equation (130) can also be expressed as a function of k=(2miEkin)1/2/. Additional temporal averaging of the carrier distribution function Eq. (130), for example over a simulation subinterval Δ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 tsim independently of the chosen initial conditions. However, the required tsim 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 tf. For a time independent scattering rate τ01, the free flight time is given by tf=τ0ln(r), where r is a random number evenly distributed between 0 and 1.149 Here, we choose τ01mjW̃ik,j(m), corresponding to an upper bound for the total outscattering rate Wik(t)=mjWik,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 Δ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 τ0Wik,j(m)(t). The too high amount of scattering assumed is corrected by introducing artificial “self-scattering,” which occurs with a probability 1τ0mjWik,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. The final wave vector magnitude k 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 Ne electrons are located in the central period, i.e., Np = 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 ΔN of electrons scattered from the central period to the left- and right-neighbouring period, respectively, over a simulation subinterval Δt, the current density can be computed as

J=ΔNeΔtS=ΔNensNpΔtNe.
(131)

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 and |jk, respectively. The dipole matrix element then becomes

dik,jk=eik|z|jk=dijS1Sexp[i(kk)r]d2r,
(132)

with

dij=ei|z|j=eψi(z)zψj(z)dz,
(133)

and for |dik,jk|2, we obtain

|dik,jk|2=4π2|dij|2S1δ(kk).
(134)

Here, we have used that |Sexp[i(kk)r]d2r|2 can be approximated by 4π2Sδ(kk) for sufficiently large in-plane cross sections S. The photon-induced transition rate given in Eq. (54) now becomes

Wik,jkopt=4π3ϵ0cn02S|dij|2δ(kk)mImLij(ωm,k),
(135)

where we sum over all relevant cavity modes with frequencies ωm and intensities Im to account for multimode lasing. Adapting Eq. (51) to the present case, the definition of the Lorentzian lineshape function becomes

Lij(ω,k)=1πγij(k)γij2(k)+[ω|ωij(k)|]2,
(136)

where ωij(k)=(EikEjk)/ denotes the resonance frequency and γij(k) is the optical linewidth of the transition, which is given by

γij(k)=12(iWik,+jWjk,),
(137)

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

Wik,jopt=πϵ0cn02|dij|2mImLij(ωm,k).
(138)

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, i.e., nis=1/S, is obtained from Eq. (52) by summing over the transitions to all available final states |jk. With Eqs. (61) and (134), we obtain

g(ω)=πωϵ0cn0SLpjωij(k)|ωij(k)||dij|2Lij(ω,k).
(139)

For Eik<Ejk, we have ωij(k)/|ωij(k)|=1, 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 of the central simulation period(s)76 

g(ω)=πωϵ0cn0SLpnjωinj(k)|ωinj(k)||dinj|2Linj(ω,k),
(140)

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 ns given by S=Ne/(nsNp), where NP 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 γij has to be taken in Eq. (136). Furthermore, ωij does not depend on the wave vector if nonparabolicity effects are neglected. Then the gain simplifies to Eq. (52), thus becoming

g=πωϵ0cn0Lpi,j|dij|2(nisnjs)Lij(ω).
(141)

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

n0c1tIm=Γmg(ωm)ImamIm.
(142)

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).

FIG. 17.

Illustration of coupled simulations including both the carrier transport and optical cavity field.

FIG. 17.

Illustration of coupled simulations including both the carrier transport and optical cavity field.

Close modal

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.

FIG. 18.

Simulated unsaturated and saturated power gain coefficient vs. frequency. The dashed line indicates the threshold gain. Reprinted with permission from Appl. Phys. Lett. 96, 011103 (2010). Copyright 2010 AIP Publishing LLC.76 

FIG. 18.

Simulated unsaturated and saturated power gain coefficient vs. frequency. The dashed line indicates the threshold gain. Reprinted with permission from Appl. Phys. Lett. 96, 011103 (2010). Copyright 2010 AIP Publishing LLC.76 

Close modal
FIG. 19.

Comparison of the measured optical power vs. lattice temperature178 to EMC simulation results; the inset contains the simulated power vs. applied bias for various lattice temperatures. Reprinted with permission from Appl. Phys. Lett. 96, 011103 (2010). Copyright 2010 AIP Publishing LLC.76 

FIG. 19.

Comparison of the measured optical power vs. lattice temperature178 to EMC simulation results; the inset contains the simulated power vs. applied bias for various lattice temperatures. Reprinted with permission from Appl. Phys. Lett. 96, 011103 (2010). Copyright 2010 AIP Publishing LLC.76 

Close modal

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 efficiency9 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.

FIG. 20.

The measured current-voltage characteristics for a high-efficiency mid-infrared QCL9 is compared to results obtained from EMC simulations with and without lasing included. Reprinted with permission from J. Appl. Phys. 110, 013108 (2011). Copyright 2011 AIP Publishing LLC.106 

FIG. 20.

The measured current-voltage characteristics for a high-efficiency mid-infrared QCL9 is compared to results obtained from EMC simulations with and without lasing included. Reprinted with permission from J. Appl. Phys. 110, 013108 (2011). Copyright 2011 AIP Publishing LLC.106 

Close modal

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 Δ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 Ω=ΔE/,125 and the corresponding tunneling time corresponds to half the duration of an oscillation cycle, τtun=π/Ω.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 ττ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 ΔE and thus long tunneling times τtun.125 Especially in terahertz QCLs, thick injection barriers are used with ΔE1 meV,125 corresponding to τ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.

FIG. 21.

Conduction band profile and probability densities for a terahertz QCL, computed based on (a) the actual potential V modeled as a periodic sequence of stages and (b) the tight-binding conduction band profile Vtb obtained by extending the barriers at the stage boundaries to confine the wavefunctions within the stage. The rectangles denote a single stage. The extended wavefunctions spanning the thick barrier and the corresponding localized wavefunctions are marked by bold lines.

FIG. 21.

Conduction band profile and probability densities for a terahertz QCL, computed based on (a) the actual potential V modeled as a periodic sequence of stages and (b) the tight-binding conduction band profile Vtb obtained by extending the barriers at the stage boundaries to confine the wavefunctions within the stage. The rectangles denote a single stage. The extended wavefunctions spanning the thick barrier and the corresponding localized wavefunctions are marked by bold lines.

Close modal

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 Vtb 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 Ωij=1i|VVtb|j.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 ψ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 ρij (see also Sec. V B), which can represent both pure and mixed states and thus is suitable to model dephasing. The diagonal terms ρii correspond to the occupation of level i, and the off-diagonal elements ρ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, ρii=nis. The time evolution of the density matrix ρij(t) is described by the von Neumann equation

idtρij=(HiρjρiHj),
(143)

where the Hamiltonian matrix elements are defined as Hij=i|Ĥ|j. For describing the carrier transport in QCLs, the Hamiltonian is divided into two parts, Ĥ=Ĥ0+Ĥ. Here, Ĥ0 describes the coherent evolution of the quantum system due to the conduction band potential, while Ĥ contains the perturbation potentials corresponding to the various scattering mechanisms and introduces dissipation to the system. Ĥ 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 form111,113

dtnis=jiτji1njsτi1nis+jiΩij(ρijρji),dtρij=iΩij(nisnjs)iρijωijγijρij,
(144)

where ωij=(EiEj)/, and τi1=jiτij1 indicates the total inverse lifetime of level i. The scattering rates τij1 have already been discussed in Secs. V A and VI E. The dephasing rate is given by125 

γij=(τi1+τj1)/2+τpure,ij1.
(145)

Here, (τi1+τj1)/2 is the lifetime broadening contribution [see also Eq. (137)], and τpure,ij1 contains the pure dephasing. The Rabi frequency Ωij is nonzero only for doublets spanning a coupling barrier; specifically, Ωii=0. In Eq. (144), resonant tunneling is assumed to be independent of the in-plane wave vector, and the rates τij1,γij are averaged over the kinetic electron distribution within the subbands.111 The stationary solution of Eq. (144) is obtained by setting dt = 0, yielding113 

ji[τji1njs+Rij(njsnis)]τi1nis=0,
(146)

with

Rij=2Ωij2γij11+ωij2γij2.
(147)

Furthermore, Eq. (38) has to be fulfilled, i.e., the total sheet density in each period is determined by the doping sheet density ns.

Equation (144) can be solved using empirical rates, or implemented in a self-consistent manner by calculating the scattering rates τij1 as discussed in Sec. VI E. The pure dephasing rate τpure,ij1 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 τ121=τ211=0 to τ21n2s+R12(n2sn1s)=0. The current density through the barrier is then given by J=eR12(n1sn2s)=eτ21n2s. With Eq. (147) and n1s+n2s=ns, we obtain10 

J=2ensΩ122γ1211+ω122γ122+4Ω122γ121τ2.
(148)

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 rather than |i. Consequently, the density matrix is given by ρik,jk, and i,j, in Eq. (143) have to be replaced by ik,jk, and k, 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 

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 mass200,201 and k·p method,202 to atomistic representations such as tight binding203,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.

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

Ĥ0=22z1m(z,E)z+2k22m(z,E)+V(z),V(z)=Vc(z)eΦ(z),
(149)

where k is the lateral electron momentum, Φ(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 k209 

GR=(E1̂Ĥ0ΣR)1,
(150a)
G<=GRΣ<GR,
(150b)
Σ<=G<D<,
(150c)
ΣR=GRDR+GRD<+G<DR.
(150d)
Here, ΣR and Σ< denote the retarded and lesser self-energies and DR and 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 Ĥ0 do not increase the numerical complexity of NEGF. The Green's functions and self-energies are functions of two spatial coordinates z,z, 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)GR(z,z,0,E)] shows width and location of resonant states in the system. This indicates that GR contains the information of resonant states and the density of states. The 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,

n(z)=dEn(z,E)=2(2π)3dEd2kG<(z,z,k,E),
(151)
J(z)=dEJ(z,E)=e(2π)3limzzdEd2k1m(z,E)×(zz)G<(z,z,k,E).
(152)

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)

Ĵ(x1,t1)=limt2t1limx2x1iv̂G<(x1,t1;x2,t2),
(153)

with the velocity operator

v̂=i[Ĥ,x̂].
(154)

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 α(z,ω) for a photon of frequency ω is a function of the permittivity ϵ(z,ω)212 

α(z,ω)=[ϵ(z,ω)ϵ0]2ωc{[ϵ(z,ω)ϵ0]+|ϵ(z,ω)ϵ0|}1/2.
(155)

Here, ϵ0 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α. The permittivity depends on the complex conductance σ(z,ω) and the materials dielectric constant ϵr(z)

ϵ(z,ω)=ϵ0ϵr(z)+iσ(z,ω)/ω.
(156)

Before lasing starts and the perturbation δV̂(ω) 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 

δG<(ω,E)=GR(E+ω)δV̂(ω)G<(E)+G<(E+ω)δV̂(ω)GR(E)+GR(E+ω)δΣR(E+ω,E)G<(E)+G<(E+ω)δΣR(E+ω,E)GR(E)+GR(E+ω)δΣ<(E+ω,E)GR(E).
(157)

Changes to the scattering self-energies δΣ due to the photonic perturbation δV̂ 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

δV̂(ω)=em(z)ωEz(ω)z,
(158)

with the photon electric field component in z direction Ez. This results in the change of the current density δJ(z,ω) in linear order of δV̂(ω)

δJ(z1,ω)=limz2z12em(z1)(z1z2)×dE2πd2k(2π)2δG<(z1,z2,k,ω,E)2e2Ez(ω)m(z1)ωdE2πd2k(2π)2G<(z1,z1,k,E).
(159)

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 δG<. The second term results from the change of the velocity and of the current operator due to the fact that δV̂(ω) does not commute with the position operator [see Eq. (154)]. The quotient of the perturbation of the current density δJ(z,ω) and the electric field Ez(ω) of the photon gives us the optical conductance

σ(z,ω)=δJ(z,ω)Ez(ω).
(160)

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

σ(z1,ω)=limz2z12e2m(z1)2(2π)3ω(z1z2)dEd2kdz3×[GR(z1,z3,k,E+ω)zG<(z,z2,k,E)|z=z3+G<(z1,z3,k,E+ω)zGA(z,z2,k,E)|z=z3]2e2dEd2k(2π)3m(z1)ωG<(z1,z1,k,E).
(161)

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 δΣ in Eq. (157) increase the numerical load significantly, since they require self-consistent iterations of δG with δΣ. Both δG and δΣ depend on the optical frequency ω in addition to the “standard” dependence of the Green's function [e.g., δΣ=δΣ(z,z,k,E,ω)]. 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 δΣ may not predict quantitative values of α(z,ω), 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

Boundary conditions in stationary NEGF always fall into two categories: boundary conditions for GR 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 GR is van-Hove singular whenever the energy 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 GR and G< are then introduced by including contact self-energies ΣconR and Σcon< in the equations for GR and 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 Σ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

f(z,k,E)iG<(z,z,k,E)/A(z,z,k,E),
(162)

with the spectral function

A(z,z,k,E)=i[GR(z,z,k,E)GR(z,z,k,E)].
(163)

In equilibrium, the function f(z,k,E) equals the Fermi distribution f(E,μ).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

f(0,k,E)=!f(nLp,k,E+neϕ),f(L,k,E)=!f(LnLp,k,Eneϕ),

field periodic boundary conditions for G< are mimicked.122 Here, Lp is the length of a QCL period, ϕ 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.

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-energies195,209,211

Σ<(z3,z4,k,E)=1(2π)3d2qdqz|CQ|2eiqz(z3z4)×[NQG<(z3,z4,|kq|,EωQ)+(1+NQ)G<(z3,z4,|kq|,E+ωQ)],
(164)
ΣR(z3,z4,k,E)=1(2π)3d2qdqz|CQ|2eiqz(z3