Effective discrete-level density matrix model for unipolar quantum optoelectronic devices

1 Introduction

Increasingly, quantum confinement in semiconductor heterostructures is exploited to develop quantum optoelectronic devices with enhanced performance and expanded functionalities. In unipolar devices, the lasing transition occurs between quantized states in the conduction band, and thus the optical properties do not depend on the semiconductor bandgap. This opens up enormous possibilities for custom-tailoring lasing wavelengths, optical nonlinearities and other active region properties by quantum engineering the confined states. Specifically, the quantum cascade laser (QCL) utilizes optical intersubband transitions in the conduction band to access a wide range of mid-infrared (MIR) and terahertz (THz) wavelengths [1], [2]. Here, a periodic multi-quantum well design is used, allowing for the generation of multiple photons by a single injected electron. Also amplifiers [3], [4], modulators [5] and detectors [6], [7] have been realized based on this principle. Generally, unipolar quantum well devices have an enormous potential for long-wavelength optoelectronic applications [8]. Furthermore, also semiconductor quantum wire structures with two-dimensional quantum confinement are attractive candidates for developing intersubband optoelectronics [9].

Recently, dynamic waveform generation with unipolar devices has become a vibrant research field, motivated by a wide range of applications in, e.g., metrology and communications. In particular, mode-locking techniques have been employed for generating short-pulse trains [10], [11] and broadband frequency combs [12], [13], [14], i.e., discrete, equally spaced spectra associated with periodic temporal waveforms. In this context, also harmonic operation in QCLs has attracted considerable interest, where the waveform period is a harmonic of the cavity roundtrip time [15], [16], [17], [18], [19], [20]. Moreover, the formation of dissipative solitons in QCLs has recently caught wide attention [21], [22], [23], [24]. For a systematic design of such waveform-generating nanostructured lasers and improved understanding of their dynamics, accurate and efficient numerical models are essential [25]. More generally, such dynamic modeling approaches are potentially relevant for high-speed systems employing unipolar quantum optoelectronic devices. To account for quantum coherence effects, these approaches are frequently based on a density matrix (DM) formalism describing the electron dynamics in the quantum active region, coupled to Maxwell’s equations capturing the optical field propagation in the cavity. Since often simulations over many hundred or thousand cavity roundtrips are required to reach steady state operation [26], the model is commonly simplified to reduce the numerical load. For example, optoelectronic devices with a waveguide cavity typically feature an invariant transverse field distribution, enabling the use of a one-dimensional optical propagation model which only depends on time t and a single spatial coordinate x [25]. Furthermore, the dependence on the electron in-plane wavevector k is typically ignored in the dynamic DM equations [25], greatly reducing the numerical load in comparison to fully k dependent models [27], [28], [29], [30]. This is justified for optical transitions between subbands with nearly parallel dispersion relationships [25], [31], as is often ideally assumed for QCLs [1], but not for interband transitions since the energy dispersions in the conduction and valence bands have opposite curvatures [25], [32]. However, also operation in unipolar quantum well and wire devices can be affected by residual nonparabolicity [33], [34], [35], [36] as well as by Bloch gain [34], [37], [38], [39], [40], both leading to gain asymmetry and linewidth enhancement.

Restricting the description of the quantum active region to two energy levels results in the semi-phenomenological Maxwell-Bloch (MB) equations, which include dissipation in terms of empirical relaxation rates [25]. Various strategies have been employed to derive effective MB equations for bipolar semiconductor lasers and amplifiers from microscopic models by adequate k averaging over the electron and hole distributions [32], [41], [42]. These models include a linewidth enhancement factor (LEF) to describe nonparabolicity effects. For unipolar devices, the so-called effective semiconductor MB equations (ESMBEs) have been derived by combining the MB equations with a phenomenological expression for an asymmetric material susceptibility [43], [44], and employed for studying the dynamic QCL operation in both ring and Fabry–Perot configurations [35], [43], [44], [45], [46]. Also the Bloch gain has been implemented in the MB framework [39]. On the other hand, fully quantitative modeling of quantum-engineered optoelectronic devices requires explicit consideration of all relevant mechanisms and quantized energy levels. This can be achieved in the framework of an advanced Maxwell-DM model, featuring a multilevel DM and a generalized system Hamiltonian, which generally includes tunneling in addition to light–matter interaction [25], [47]. Dissipation is here described using the Lindblad formalism [25], [48]. The Lindblad-type relaxation terms and Hamiltonian matrix elements can be extracted from carrier transport simulations or microscopic descriptions, resulting in a self-consistent device model [26], [31]. This approach has been employed for quantitative simulations of various advanced THz and MIR QCL devices in Fabry–Perot and ring configurations, yielding excellent agreement with experiment and providing insights into device operation. Examples include the modeling of soliton generation [22], short-pulse mode-locked operation [11], [47], and fundamental [26], [31], harmonic [20], [49] as well as difference-frequency comb [50] generation. The multilevel DM naturally includes gain asymmetry due to multiple optical transitions, which can have a significant influence on the optical dynamics [17]. However, contributions of nonparabolicity and Bloch gain have to date not been considered in Maxwell-DM approaches beyond the two-level approximation. In the present work, these effects are systematically incorporated by suitable k averaging of the microscopic DM equations. The resulting correction factors, such as effective transition frequencies and LEF-related quantities, are in our approach not treated as fitting parameters, but can be extracted from carrier transport simulations together with the other required parameters. Thus, the resulting effective Maxwell-DM equations preserve the self-consistent nature of the simulation model.

2 Microscopic model

For interband transitions, the derivation of effective two-level models by suitable wavevector summation has been addressed in previous work [32], [41], [42]. Here, we focus on unipolar devices. As illustrated in Figure 1, these utilize optical transitions between quantized energy levels n, each consisting of a quasi-continuum of states n , k with eigenenergies E _n,k . The associated transition frequencies are given by ω m n , k = E m , k − E n , k / ℏ , with the reduced Planck constant ℏ. Specifically for quantum well structures, quantum confinement in growth direction z results in the formation of quantized states n, and the free in-plane carrier motion is described by the two-dimensional in-plane wavevector k. Nonparabolicity can be included by allowing for an energy dependent effective mass when solving the Schrödinger equation. This yields for each subband n the corresponding wavefunction ψ n z and the electron dispersion relation

with the subband effective mass m n * and E _n = E _n,k=0 [51], [52]. Here, we assume decoupling between the in-plane motion and the confinement direction implying k independent wavefunctions ψ n z , which is a good approximation for not too narrow finite quantum wells [53]. Within this model, the effect of nonparabolicity is accounted for by the different value of m n * for each subband. We note however that the treatment of nonparabolicity in our effective DM equations derived in Section 3.2 is not restricted to dispersion relations of the form Eq. (1).

Figure 1:

Schematic representation of level schemes for unipolar devices. Additionally, parabolic electron dispersion relations as given in Eq. (1) are illustratively sketched. The upper and lower levels of the optical transition are indicated by blue and red colors, respectively. For QCLs, a periodic repetition of identical stacks (marked by rectangles) is used.

The DM elements are given by m , k ρ ̂ n , k = ρ m n , k , and reduction to an effective discrete-level description using DM elements of the form ρ _mn requires suitable k averaging. Also semiconductor quantum wire structures with two-dimensional quantum confinement are attractive candidates for developing intersubband devices [9]. They can be described analogously; in this case, the free carrier motion in the remaining direction is characterized by a one-dimensional wavevector [25]. Besides, our approach is not restricted to the parabolic dispersion relations given in Eq. (1).

The diagonal DM elements ρ _nn,k = ρ _n,k can be written as

where the distribution function f _n,k gives the electron occupation probability of a state n , k . The scaling factor S may be chosen such that ρ _n = S∑ _k f _n,k corresponds, e.g., to the carrier number density in level n, or that the normalization condition ∑ρ _n = 1 is fulfilled. The off-diagonal elements ρ _ij,k contain the coherence between states i , k and j , k . The DM evolution equation is given by

with the collision term ∂ t ρ ̂ col . DM elements between different wavevectors need not be considered due to the k conservation of optical transitions. This also applies for first-order tunneling processes, which can straightforwardly be included in the Hamiltonian by employing adequately localized basis states [54], [55]. For an orthogonal basis set, we obtain from Eq. (3)

The Hamiltonian in Eq. (4) can be represented as H ̂ = H ̂ 0 + H ̂ I , where H ̂ 0 is the Hamiltonian of the unperturbed system, with H _0,nn,k = E _n,k . Off-diagonal elements can, e.g., arise from the inclusion of resonant tunneling. For example, in QCL designs electron transport across thick injection or extraction barriers is mediated by tunneling between closely aligned states. Using a localized basis, such as tight-binding states m , k and n , k located at the left and right of the barrier [55], [56] or EZ states [57], the corresponding DM elements can for small ω m n , k be written as H _0,mn = H _0,nm = ℏΩ _mn . Assuming k independent wavefunctions as discussed above, also the coupling energy ℏΩ _mn = ℏΩ _nm does not depend on k [55], [56]. For an optical transition between two states m , k and n , k , light–matter interaction can in dipole approximation be described by the corresponding matrix elements of the interaction Hamiltonian, H _I,mn = H _I,nm = −Ed _mn , with the optical field E t . Here, d _mn = d _nm represents the dipole matrix element, which is again k independent under above assumptions. The model for the collision term in Eq. (4) used here is discussed in Appendix A. We note that within above framework, second-order effects connecting states with different wavevectors, such as Bloch gain and second-order tunneling, are not yet included. Various approaches have been discussed in literature to consider these contributions in effective discrete-level DM models [39], [58], [59]. In Section 3.2.1, we give a detailed discussion on the implementation of Bloch gain.

In the following, we restrict our discussions to a field E with moderate bandwidth and in close resonance with the optical transition(s), i.e., ω c ≈ ω m n , k where ω _c denotes the optical carrier frequency. Furthermore assuming non-excessive field strengths as is typically justified in optoelectronics, the widely used rotating wave approximation (RWA) can be invoked to increase the numerical efficiency of the model [25]. Here, the fast oscillations of E and the related DM elements ρ _mn,k around ω _c can be removed by representing these quantities in terms of the slowly varying envelope functions ɛ _mn and η _mn,k ,

More specifically, the field envelope ɛ _mn = ɛ _nm is here expressed in terms of the corresponding (instantaneous) Rabi frequency. The asterisk denotes the complex conjugate, and sgn represents the sign function. The evolution equations for the DM elements in RWA are then obtained in the usual manner by substituting Eqs. (5a) and (5b) into (4) and discarding the rapidly oscillating terms (see Appendix B). Since a coarser spatiotemporal grid can be used to resolve the dynamics of the envelope functions, the computational load gets significantly reduced as compared to full-wave simulations.

3 Effective discrete-level model

For the DM-based dynamic modeling of semiconductor lasers and other optoelectronic devices, typically a two- or multilevel model featuring discrete energy levels is used, where the wavevector dependence of the states is not explicitly taken into account [25]. Besides the considerable decrease of numerical complexity as compared to fully microscopic models [25], discrete-level approaches facilitate the development of compact and intuitive descriptions of the laser dynamics [45], [60], [61], [62]. However, the wavevector dependence may leave a direct imprint on the DM dynamics beyond microscopic interactions, e.g., in form of an asymmetric susceptibility and the closely related linewidth enhancement resulting from the nonparabolicity effect or Bloch gain [32], [35], [39], [43]. Thus, rather than simply ignoring the wavevector dependence of the microscopic states, a systematic removal of this quantity from the model by adequate k summation is more appropriate.

The transition from the microscopic, k-resolved description to an effective model is achieved by defining effective DM elements obtained via k summation,

where the diagonal DM elements, ρ _nn ≕ ρ _n , are related to the total population of level n, and the elements η _mn to the polarization of the optical transition m → n. For a stationary optical field with frequency ω, η _mn is directly proportional to χɛ _mn with the complex susceptibility χ ω . In Figure 2, χ is schematically illustrated for the case of parallel subbands and for nonparabolicity, resulting from different effective masses of the upper and lower subband. Here, population inversion is assumed. For the complex field convention introduced in Eq. (5b), I χ is proportional to the loss coefficient. For parallel subbands, both the harmonic gain and the real part of the Bloch susceptibility assume the typical Lorentzian shape. An asymmetric susceptibility is obtained for nonparabolicity, or also for parallel subbands if both the harmonic and Bloch contributions are present.

Figure 2:

Schematic representation of harmonic and Bloch contribution to the susceptibility χ for parabolic subbands and for nonparabolicity. Here, ω ₂₁ and γ ₂₁ are the resonance frequency at k = 0 and the dephasing rate.

4 Examples

Nonparabolicity is usually much more pronounced in mid-infrared (MIR) than in terahertz QCLs, since the larger energy spacing between the upper and lower laser levels tends to enhance the difference between the effective masses. Furthermore, the nonparabolicity effect increases with electron temperature since higher k states get occupied. Thus, in the following we focus on high-temperature MIR QCL structures.

4.1 Analytical effective parameter model

In order to validate the analytical effective parameter model introduced in Section 3.2.3, we choose m u * = 1.2 m ℓ * = 0.06 m e , ρ _ℓ = ρ _u /3, T _ℓ = 1.5 T _u = 900 K, γ _uℓ = 10 ps⁻¹ corresponding to a Lorentzian gain bandwidth of γ _uℓ /π = 3.2 THz, and γ _u = γ _ℓ = γ _uℓ /2. These are realistic values for MIR QCLs and give rise to a pronounced nonparabolicity. From Eqs. (18), (29) and (31), the effective parameter values c u ℓ , u np = 0.886 + 0.228 i , c u ℓ , ℓ np = 0.796 + 0.248 i , c u ℓ , u b = − 0.013 + 0.054 i , c u ℓ , ℓ b = 0.012 − 0.035 i , and Γ _uℓ = 14.2 − 5.2i ps⁻¹ are obtained. For validating the effective model, we compare the frequency dependent susceptibility χ ∝ η _uℓ /ɛ _uℓ computed from Eq. (21) with the result of the fully k dependent calculation, obtained by solving Eq. (11) in analogy to Eq. (21) and employing Eq. (6). In Figure 3(a), the obtained susceptibility is shown for the effective and fully k dependent model as well as for the conventional discrete-level DM equations, obtained by setting Γ _uℓ = s _uℓ and c _uℓ,u = c _uℓ,ℓ = 1 in Eq. (21). Figure 3(b) displays the harmonic and Bloch contributions to χ. Overall, we find good agreement between the full and the effective model. As expected, the conventional discrete-level DM equations do not capture the asymmetry and broadening of χ and α caused by nonparabolicity and the Bloch contribution.

Figure 3:

Calculated (a) susceptibility χ and (b) harmonic and Bloch contributions to χ as a function of the normalized frequency detuning Δ/γ _uℓ for a two-level system. The results from the analytical effective parameter model, introduced in Section 3.2.3, are compared to calculations based on the conventional discrete-level and the fully k resolved DM model.

4.2 Multilevel effective parameter model

As a test case for the general effective multilevel DM model of Section 3.3, we choose a diagonal bound-to-continuum room temperature QCL design emitting at 8.5 μm [66], which has been widely used as a reference structure for validating modeling approaches [65], [66], [67]. In Figure 4(a), the energy levels of a representative stage, which have been computed with a Schrödinger–Poisson solver, are displayed. Furthermore, DM-Monte Carlo carrier transport simulations have been performed [56], [65]. The simulated energy dependent distribution functions and dephasing rates are shown in Figure 4(b) and (c), respectively. The electron dispersion relation is here modeled using Eq. (1). The effective masses of the upper and lower laser levels u and ℓ are 0.0604 and 0.0547, giving rise to a pronounced nonparabolicity. Additionally, Bloch gain between the laser levels and the coherent coupling of the upper laser level to the tunneling injector t contribute to the asymmetry. For this subset of coupled subbands, Eq. (33) becomes

(38) ∂ t η u ℓ , k η t ℓ , k ρ t u , k = − s u ℓ , k i Ω t u 0 i Ω t u s t ℓ , k i 2 ε u ℓ 0 i 2 ε u ℓ * s t u , k η u ℓ , k η t ℓ , k ρ t u , k + i 2 ε u ℓ ρ ℓ , k − ρ u , k 0 i Ω t u ρ t , k − ρ u , k ,

with ρ _i,k = Sf _i,k where the scaling factor S introduced in Eq. (2) can be freely chosen. Taking advantage of the in-plane isotropy, we represent the k dependence in terms of the energy variable w introduced in Eq. (28). The ρ i w and s i j w are provided by the carrier transport simulations at the operating point as shown in Figure 4; furthermore, ℏΩ _tu = 3.2 meV is obtained. Assuming operation close to threshold, an arbitrary small value for ɛ _uℓ can be assumed. Setting ∂ _t = 0, a linear equation system is obtained for σ i j o w = η u ℓ w , η t ℓ w and ρ t u w . Plugging the results in Eq. (35) and replacing the summation over k by integration over w, the correction coefficients Γ _mn,ij are obtained. The reduced effective DM equations are then given by

(39) ∂ t η u ℓ η t ℓ ρ t u = − Γ u ℓ i c u ℓ , t ℓ Ω t u 0 i c t ℓ , u ℓ Ω t u Γ t ℓ i 2 c t ℓ , t u ε u ℓ 0 i 2 c t u , t ℓ ε u ℓ * Γ t u ⋅ η u ℓ η t ℓ ρ t u + i 2 ε u ℓ c u ℓ , ℓ ℓ ρ ℓ − c u ℓ , u u ρ u 0 i Ω t u c t u , t t ρ t − c t u , u u ρ u ,

where c m n , i j = Γ m n Γ m n , i j − 1 , and Eq. (37) has been used to include Bloch gain.

Figure 4:

Carrier transport simulation results for the investigated QCL. (a) Conduction band profile with energy levels and probability densities. (b) Electron distribution functions f i w for the upper laser level (u), lower laser level (ℓ), and tunneling injector (t) as a function of the energy w, defined in Eq. (28). (c) Dephasing rates γ i w and γ i j w .

For validating Eq. (39), we again compare the frequency dependent susceptibility χ ∝ η _uℓ /ɛ _uℓ . Similarly as for Eq. (21), the frequency dependent η _uℓ is obtained by inserting a frequency-shifted field ε u ℓ → ε u ℓ Δ exp − i Δ t and corresponding optical DM elements η i j → η i j Δ exp − i Δ t into Eq. (39), where Δ = ω − ω _c. Cancelling the factor exp − i Δ t and setting ∂ _t = 0 yields a linear equation system for the stationary solution, which is solved in dependence of Δ to obtain η u ℓ Δ . The exact result is obtained by computing η u ℓ , k Δ in an analogous manner from Eq. (38), and performing the k summation according to Eq. (6). As can be seen from Figure 5(a), the susceptibility obtained with the effective discrete-level DM model agrees well with exact result of the fully k dependent DM simulation, while the conventional discrete-level DM equations, obtained by setting Γ _mn = s _mn and c _mn,ij = 1 in Eq. (39), do not provide a good fit. The results shown in Figure 5(b) are for the same device, but the carrier transport simulations have been performed under lasing conditions [65], resulting in gain saturation. Again, the effective DM model yields good agreement with the exact results.

Figure 5:

Susceptibility χ for a three-level system, consisting of the laser levels and a tunneling injector, as a function of the frequency detuning Δ / 2 π for the (a) unsaturated and (b) saturated case. The results from the generalized effective parameter model, introduced in Section 3.3, are compared to calculations based on the conventional discrete-level and the fully k resolved DM model.

4.3 Dynamic simulations

Finally, we present simulations based on the effective Maxwell-DM approach to assess the numerical performance of the model, and to investigate the influence of nonparabolicity and Bloch gain on the QCL dynamics. For the dynamic simulations, Eqs. (C.3)–(C.9) in Appendix C are solved on a spatiotemporal grid, using an explicit 3rd order Adams-Bashforth method for Eqs. (C.3)–(C.7) and a finite difference scheme for Eq. (C.8) [25]. To obtain realistic results, SHB and group velocity dispersion are included in the model. Furthermore, spontaneous emission noise is considered in Eq. (C.8) to account for the associated field fluctuations and to emulate the buildup of lasing. As an exemplary structure, we choose a Fabry–Perot cavity with a vertical two-phonon resonance active region, featuring room temperature operation at around 9 μm [68], [69]. This design has for example been used for investigating the formation of dense and harmonic multimode spectra under different driving conditions [70]. Similarly as in Figure 4(a), we model injection into the upper laser level by tunneling though the thick injection barrier. Thus, the coherent coupling between the injection, upper and lower laser levels can again be described by Eq. (39). As outlined in Section 4.2, the Hamiltonian matrix elements, scattering/dephasing rates and effective parameters are extracted from carrier transport simulations. In Figure 6(a), the computed susceptibility at lasing threshold is shown as a function of frequency for the same models as in Figure 5, again yielding excellent agreement between the effective and full DM approach. In addition, the Bloch and harmonic contributions are displayed in Figure 6(b) for the effective DM model. Both the nonparabolicity and the Bloch gain contribute to the gain asymmetry, resulting in a noticeable shift of the gain peak to lower frequencies.

Figure 6:

Simulation results for QCL multimode operation: (a) active region susceptibility χ at threshold, calculated with different models; (b) Bloch and harmonic contributions to χ according to the effective parameter model; (c) and (d) multimode spectra obtained with the (c) effective and (d) conventional Maxwell-DM approach.

In the following, we focus on dense multimode operation, since the emergence of harmonic spectra in free-running lasers is quite elusive, critically depending on the drive history, sample used and other factors [19], [49], [70]. Exemplarily, we investigate the effect of nonparabolicity, since the influence of Bloch gain on the QCL dynamics has already been studied in detail for a similar active region design [39]. In Figure 6(c) and (d), simulation results of the effective and conventional Maxwell-DM model are shown for a moderate two-facet output power of ∼ 50 m W . Both approaches yield a dense multimode spectrum already slightly above threshold, as also observed in experiment [70]. Although multimode operation in Fabry–Perot cavities is largely governed by SHB [71], [72], the spectrum is clearly broader for the effective model, featuring a 20-dB bandwidth of 0.71 THz (i.e., 2.1 % relative bandwidth) versus 0.54 THz (1.6 %) for the conventional approach. This illustrates the contribution of nonparabolicity-induced linewidth enhancement to multimode formation. In addition, the spectrum obtained with the effective model is downshifted in frequency and thus agrees somewhat better with the experimentally observed wavelength range [70], which however has been found to depend significantly on the growth process and facility [73]. Simulations at higher output powers likewise yield broader and frequency-downshifted spectra for the effective model, providing a better overall match to experiment as expected. The numerical stability and efficiency of the effective Maxwell-DM approach has been further validated by applying it to other test structures. Since the effective parameters can directly be extracted from the carrier transport simulations and the effective Maxwell-DM equations have the same complexity as the conventional model, the computational cost is comparable for both approaches. Thus, the effective Maxwell-DM model is well-suited for dynamic QCL modeling, and specifically for the investigation of operating regimes where gain asymmetry plays a pronounced role, such as comb and soliton formation in ring cavities [21], [22], [23], [24], [39], [43], [45], [74] and harmonic operation [15], [16], [17], [18], [19], [20], [70].

5 Conclusions

An effective DM model has been derived for unipolar quantum optoelectronic devices by adequate summation over the electron wavevector, which characterizes the free carrier motion in the directions without quantum confinement. The resulting effective discrete-level DM equations differ from models for true discrete-level quantum systems, such as quantum dots, by containing additional effective parameters. This extended description includes gain asymmetry and linewidth enhancement by considering effects such as nonparabolicity and Bloch gain. Here, the effective parameters are extracted from carrier transport simulations, providing a self-consistent model without phenomenological parameters. Good agreement with fully wavevector dependent simulations is found. By coupling the DM description to optical propagation equations, an effective Maxwell-DM model is obtained for the combined optical and electronic device dynamics. The approach is validated by exemplary QCL simulations, achieving numerical performance comparable to the conventional discrete-level model while offering greatly improved accuracy and versatility. Thus, the effective Maxwell-DM equations are well-suited for the theoretical investigation of dynamic operating regimes, such as comb generation in ring cavities or the formation of solitons and harmonic states. The predictive power of the model may be further enhanced by taking into account the contributions of non-resonant optical transitions to linewidth enhancement. Perspectively, an adaption of the presented approach to bipolar quantum optoelectronic devices would be highly attractive. In this context, interband cascade lasers (ICLs) [75] are of particular interest, since they have recently shown great potential for the generation of dynamic waveforms in the mid-infrared regime, such as short pulses [76], broadband frequency combs [77], [78] and harmonic comb states [79]. Suitable approaches for microscopic carrier transport simulations, required as input for the self-consistent dynamic device model introduced in this paper, are meanwhile available for ICLs [80]. Generally, for bipolar optoelectronic devices a main challenge is that computing the effective parameter integrals may involve divergence problems [32], which must be adequately handled.