A density matrix model of transport and radiation in quantum cascade lasers

A transport model for quantum cascade lasers based on density matrix formalism that incorporates the laser optical field is confronted with experiment. For a typical mid-infrared laser, very good agreement is found for both the current–voltage and current–optical power characteristics. Forcing thermal distribution with a unique temperature in all subbands was found to lead to an overestimate of electron heating in the injector. The model can then be used further to optimize and design new structures.


Introduction
Since their invention [1], quantum cascade lasers have been modelled in various ways, from rate equations [2,3], density matrix [4]- [7] and Monte-Carlo [8]- [11] simulations to pure quantum mechanical models [12]- [15] that account fully for the in-plane dynamics. The simplest models, based on a few fundamental parameters [16], were able to predict the threshold current density accurately but not the light-current-voltage curves. It turns out that innovation in design requires predictive simulation tools. Between a simple scattering model and full quantum mechanical simulations, there is room to develop effective models that enable both quantum effects, with sequential resonant tunnelling [6], and account for the scattering of electrons through rate equations. These models have the advantage of being light enough to allow computation of the laser optical field and its effect on the transport.

Choice of an eigenbasis
In quantum cascade lasers, the choice of a basis of eigenstates on which transport can be computed is crucial. Indeed, the Fermi golden rule will only hold if the energy spacing between eigenstates E is much larger than broadening γ of the states themselves. While traditionally the eigenstates were computed on a one period length, in our model in order to keep a large enough E, we use a basis that spans a portion of each period of the structure. As shown in figure 1, the period is divided into two sub-periods: the active wells region that includes the laser doublet and the two depletion states, and the injector region designed to ensure the extraction and the relaxation of the electrons to the next stage. At the splitting barriers, the transport is  Figure 1. Two active periods shown at the injection resonance (44 kV cm −1 ). The coupling between the sub-periods, achieved by resonant tunnelling (coupling energy matrix for injection and˜ for extraction), is shown through injection/extraction barriers. The rate equations from the upper laser state are represented in the active wells of the second period. Light pink layers are doped. modelled by sequential resonant tunnelling. A simplified model that explicitly incorporates the resonant nature of the transport has been successfully applied in terahertz quantum cascade lasers [5]. Another important aspect is the second-order contributions to tunnelling that have been outlined for the gain [17,18] and current [6] between subbands.

Structure of the model
We have developed and numerically implemented a general density matrix model that allows splittings at any barriers in the structure. In each sub-period, the Fermi golden rule is applied to compute scattering between subbands. We do not solve the problem completely in k-space, but we average the scattering time with a Fermi-Dirac distribution. We can therefore reduce the electron dynamics to a set of rate equations by computing the averaged intersubband broadening [19] inter with interface roughness [20], LO-phonons, alloy disorder and ionized impurities (dopants).

Resonant-tunnelling contribution to rate equations
The density matrix problem [5,6,10] can be summarized into two main equations, firstly for the populations ρ ii = n i and secondly for the coherences ρ i j : where W ik is the averaged transition rate from state k to state i,h ik is the coupling energy between state i and state k (it is nonzero only when the doublet spans a coupling barrier),h ik is the detuning between these states and τ ||,i j is the dephasing time. The effective parameters σ i j account for electrons tunnelling at a constant energy rather than at a constant wavevector [6], which effectively implements the second-order scattering in the current calculation. These equations can be solved in steady-stateṅ i =ρ i j = 0, and after some algebraic manipulations the equation of motion can be cast into the form The resonant-tunnelling contribution to rate equations reads The dephasing time that relates the loss of the phase correlation between two states is computed using equation (11) in [10]. It has two contributions. The first one comes from intersubband transitions, when an electron is scattered from a state of the resonant doublet to another subband. The second one is a pure phase contribution that we identify as intrasubband broadening using equation (4) in [19]. We have evaluated intra at the characteristic energy of the electron distribution k B T, where k B is the Boltzmann constant and T the electronic temperature.

Current
The current in the model is computed at the coupling barrier. It has two contributions. Firstly, the coherent current arises from resonant tunnelling between the states across the barrier. Secondly, direct scattering between these states achieves an incoherent current. As we work in a basis localized at least to one period, the incoherent current is a marginal contribution owing to the weak overlap of the subband wavefunctions. For the structure investigated here, the incoherent current is about 1-3% of the total current.

Computation of the population and electronic temperature
The transport populations and the corresponding shape of the potential are computed through a few self-consistent iterations starting from an initial subband configuration where electrons are uniformly distributed. The electronic temperature is computed by solving a global kinetic energy balance equation [21]. We implicitly assumed that the electron-electron interaction is sufficiently strong to thermalize all the subbands to the same temperature as some authors did [22].

Introduction of the optical field
A photon flux density S at a fixed energyhω 0 is then included in the model. To simplify, we consider two levels i and j in a sub-period; their equations of motion are modified by simulated emission and absorption: where g i j is the total net gain between levels i and j. If we assume that g i j is a linear form in the population distributions f i (k) of the subbands-this is the case for the first-order (Lorentzian) and the second-order [17] gain formula-with n i = D k k f i k , where D k is the density of states, we can write where the notation g i j [·; ·] denotes the functional dependence of the gain in the population distributions. This property allows us to define gain cross sections g c,i i j and g c, j i j as when n i = 0 and n j = 0, and g c,i i j = 0, g c, j i j = 0 otherwise. The gain between upper subband i and lower subband j can therefore be rewritten as This linear form is particularly suited for rate equations and we can write If we make the assumption that g c,i i j , g c, j i j , g c,i ji and g c, j ji are computed for initial populations n 0 i and n 0 j , we can estimate the deviation of the initial populations, owing to photon flux density, by solving these rate equations. In particular, if this system is embedded into a self-consistent routine on the populations, it will converge to the exact solution. The main advantage of using this technique is to avoid an inhomogeneous term in the general rate equations.
Practically, the gain curve is first computed with populations found with non-radiative processes only. If we decide to let the laser work at the peak gain g p , given optical losses α for the empty waveguide and for the mirrors, we obtain a threshold condition g p − α = 0, where is the modal overlap. We then match this condition by varying S.

Comparison with experiment
In the present work, we apply the model to a reference two-phonon quantum cascade laser [23]; lattice heating can be neglected as all measurements were taken in pulsed mode. Concerning the simulation, we have taken typical parameters for the various scattering mechanisms; the interface roughness is modelled [19] with a value of = 90 Å for the correlation length of the steps, = 1.2 Å for the step height and [20] κ = 15 Å for the correlation length between interfaces. These values are not empirical but were deduced from measurements in [20,23]. They weakly influence the current-voltage characteristic, but can be crucial in the threshold current determination as interface roughness is the main broadening mechanism in mid-infrared quantum cascade lasers. The LO-phonon energy ishω LO = 32 meV. The computation time for the current and the photon flux density is about 150 s per point (at a given electric field) on In order to validate our model, we simulated the temperature dependence of the threshold current. The threshold current follows I th = I 0 exp(T /T 0 ). As a first step, we simulated the current-voltage curve using the global kinetic balance model. The result is shown in figure 2; as in [6], excellent agreement exists between experiment and theory. But the simulated threshold current is slightly higher than the experimental value. The optical power is also lower. The discrepancy increases with temperature.

Failure of the kinetic balance model
The global kinetic balance model is too pessimistic because it overestimates the electronic temperature of most subbands: at the injection resonance (J = 5.2 kA cm −2 ) with a lattice temperature of 300 K, the electronic temperature is found near 600 K. This can be explained by looking closer at the laser transition. An electron injected in the upper laser state with no excess kinetic energy will perform a non-radiative transition to the lower laser state by emitting an optical phonon (inelastic) or by elastic scattering. In any case, it will end up in the lower subband with a high in-plane momentum. The model will then try to represent these hot electrons with a thermal distribution, forcing a strong heating of all the subbands. The overestimation of the electronic temperature of the injector subbands causes backfilling in the active wells while it enables optical phonon absorption from the upper laser state to excited states. Both processes reduce the population inversion and therefore explain the increase in the threshold current.

Equilibrium between the lattice and the electrons
In figure 3, we show electron distribution versus energy at the injection resonance for two opposite cases: when the global kinetic balance is solved, and when the electronic temperature is forced to the lattice temperature. The latter case considers that the electron-lattice interaction is strong enough to thermalize the electrons-these ones essentially reside in the injector region-to the lattice temperature. This reproduces plausibly the situation of the structure shown here. In fact, we can estimate the average timeτ spent by an electron in the active wells region. For an applied electric field of 34 kV cm −1 in the middle of the dynamic range, we compute the sum of the carrier density in the subbands of the active wells region n a ≈ 0.26 × 10 11 cm −2 and the corresponding current density J ≈ 2.6 kA cm −2 . Using J = q 0 n a /2τ in [24], we haveτ ≈ 1.6 ps. During this time, the electron can emit approximatively eight optical phonons [25,26], using 200 fs as the emission rate of bulk optical phonons. It allows the electron to lose ≈256 meV of kinetic energy, covering largely the energy gap formed by the optical transition and the phonon resonances designed between the lower states of the active wells. Therefore, most electrons have no excess kinetic energy when extracted to the injector region. In the injector itself, the transport is achieved by direct scattering between a dense collection of states that allows an efficient thermalization. Therefore, we kept the subbands to the lattice temperature as shown in figure 2, at 300 K. Very good agreement is found. It is important to note that no fit parameters were used.

T 0 characteristics
Now we test the model on a larger range of temperatures by computing the T 0 parameter. The latter is robust and will show us if we have chosen the right thermal model. In figure 4(a), threshold current was plotted versus lattice temperature and then fitted with the I 0 e T /T 0 function. The value extracted from measurements is T 0 = 174 K. The simulated value with the kinetic balance model is T 0 = 111 K while the result for equal lattice and electron temperatures is T 0 = 155 K, which is in better agreement with the experiment. We expect a different situation in strained compensated structures where the electronic temperature was measured [27] clearly above the lattice temperature. The simplification made here does not hold anymore as the energy gap between the upper laser state and the injector ground state is much larger. In figure 4(b), we show the light-current curves for the lowest and the highest temperatures. The agreement is excellent in terms of threshold current and slope efficiency.

Comparison with luminescence measurements
We have successfully confronted the model with luminescence measurements at 300 K. It is important to mention that the discrepancy between the kinetic balance model and the equilibrium model does not arise from a broadening of the optical transition. We have computed the linewidth without the simplifications made in [23]: we included all scattering channels and kept the correlation terms in intrasubband calculations.

Injection efficiency
With the validation of the model, quantities that are not accessible directly by the experiment, such as injection efficiency, can be predicted. The latter can be readily computed as η inj = J upper /J total , where J upper is the current flowing into the upper laser state from the injector region and J total is the net current through the injection barrier. A value of η inj ≈ 0.86 is found.

Conclusion
We have presented a model based on density matrix formalism that enables the simulation of light-current-voltage characteristics in mid-infrared quantum cascade lasers. An important issue was the computation of the T 0 parameter. It has been found that the validity of the T 0 curve is intimately linked with the thermal model used for electrons. Our first model included subbands at the same temperature, by assuming that electron-electron interaction is strong enough to provide a thermal distribution in each subband and ensure a unique electronic temperature. Such a model was a failure, at least for the considered structure. Then we turned to an opposite model in which the subbands are in equilibrium with the lattice, yielding much better results. Anyway, to go past this model, the transport needs to be solved in k-space. This is moreover required for treating the case of local population inversion [28].
One advantage of the model is its numerical lightness, enabling automated optimization.