Numerical simulation of GaAs-based mid-infrared one-phonon resonance quantum cascade laser

Non-equilibrium Green’s function method is used to analyze electronic transport in a mid-infrared quantum cascade laser on microscopic level. Basing on the excellent agreement found between calculated and experimental data, the conclusions are derived that the carrier distribution in the lower laser subband is non-thermal, and the carriers are extracted from active region both in cold and hot state. An estimate s3 1⁄4 0:66 ps for upper lifetime at 77 K was found which considerably differs from the value 1.4 ps evaluated from the form factors.


Introduction
In-depth understanding of quantum cascade laser (QCL) operation needs advanced description that takes both scattering and quantum coherence into account. This can be achieved with density matrix (DM) or non-equilibrium Green's function (NEGF) methods (Jirauschek et al. 2014) which, however, are highly demanding, both conceptually and computationally, so their use requires simplifications. Consequently, the existing This article is part of the Topical Collection on Numerical Simulation of Optoelectronic Devices 2016.
implementations of NEGF and DM use the basis being cut to several quantum states per one QCL period (Lee et al. 2002), do not fully resolve for the in-plane momentum (Lee et al. 2002;Terazzi et al. 2010), and limit the analysis to at most three device's modules (Lee et al. 2002;Terazzi et al. 2010;Kubis et al. 2009;Hałdaś et al. 2011). They use effective mass Hamiltonians limited to at most two bands (Franckié et al. 2015), so the effect of bands mixing is also simplified. In spite of this, the results obtained with these methods agree well with the experimental data while only few parameters in the theoretical model are adjustable. This almost "fit-free" agreement validates these methods and gives the calculated characteristics, including those which are experimentally inaccessible, the status of reliable estimates of the real ones. In many cases, having such a response is crucial for the verification of the designing assumptions and better understanding of device operation. Obviously, the quality of the mentioned estimates depends on the number of the simplifications used in the model. While most of the approximations listed above can be easily released for THz devices, only few works report on the modeling of mid-IR QCLs with DM/NEGF methods that are fully energy and in-plane-momentum resolved, account for band mixing, and use numerous basis states (Hałdaś et al. 2011;Kolek et al. 2012). This paper adds to these by presenting NEGF-based analysis of the electronic transport and optical gain of mid-IR QCL emitting at ffi9 lm that uses a single-phonon resonance for lower subband depopulation.

The model
Quantum cascade lasers are purely unipolar, n-type devices, so the model uses single-band effective mass Hamiltonian. Mixing with the valence bands is taken into account, making electron effective mass, m, energy (E) dependent. Namely, mðE; zÞ ¼ m Ã ðzÞ 1 þ ½E À E c ðzÞ=E g ðzÞ È É , where E c ðzÞ and E g ðzÞ are the conduction band edge and band gap, respectively; and z is the growth direction. The in-plane dynamics is included by kinetic energy terms that assume isotropic effective mass. Such a choice preserves the inplane non-parabolicity, comparable to the results predicted by 8-band k·p method (Faist 2013). The full non-interacting Hamiltonian reads where the potential energy term V(z) comprises the conduction band edge offsets E c ðzÞ, the external bias U, and the Hartree term calculated self-consistently by the solution of the Poisson equation. Calculations are made in the real space. The Hamiltonian equation (1) is discretized and the grid points define its base vectors. This basis is quite large as the mesh of 76 grid points is used to map the QCL structure. The formulations for scattering selfenergies in this basis were provided by Lake et al. (1997) for the alloy disorder, LOphonon, interface roughness, and ionized impurity scatterings, and by Kubis et al. (2009) for the acoustic phonons scattering in energy-averaged approximation. Unlike others, the NEGF method does not a-priori assume that carriers within the QCL subbands are thermalized. Just the opposite, the intrasubband carrier distributions are calculated selfconsistently. Accordingly, there is no need to introduce some electronic temperature, and the only temperature in the model is the ambient (lattice) temperature, T. The mentioned size of the basis refers to the central period of QCL device. The remaining parts are imitated by the left and right hand contact self-energies R L;R . They were calculated by solving the Dyson equation for the QCL periods connected to the central one on its left (L) or right (R) hand side. The boundary conditions used in the calculations of RL; R assume that L/R period is connected to the semi-infinite, homogeneous leads on its left/right side and the hard wall on its right/left side. Scattering within the L/R period was included in the simplified manner by introducing dephasing, momentum relaxing (diagonal) self-energies of Golizadeh-Datta type (Golizadeh-Mojarad and Datta 2007; Kubis et al. 2009).
In the real space basis representation, the Green's functions G R ; G \ are the fourparameter functions of positions z; z 0 , the in-plane momentum modulus k, and the energy E. They can be found by an iterative solution of the Dyson, Keldysh, and Poisson equations. Then, the energy, momentum, and position-resolved densities of states (DOS), electrons (DOE), and current (J) can be found through well-known relations connecting these quantities with G R;\ ðE; k; z; z 0 Þ (Jirauschek et al. 2014;Kubis et al. 2009;Lake et al. 1997).

Results
Simulations were made for the device designed by Page et al. (2001): The GaAs/AlGaAs layer sequence 4.6, 1.9, 1.1, 5.4, 1.1, 4.8, 2.8, 3.4, 1.7, 3.0, 1.8, 2.8, 2.0, 3.0, 2.6, 3.0 nm defines one QCL period (barriers are in bold, underlined layers are n-doped). Most of the parameters in the formulations of scattering self-energies are material-dependent. They have been gathered in Table 1. In fact, the only parameters that are unknown and can be treated as adjustable are interface roughness rms height D and correlation radius K. The adjustment can be done by fitting calculated and experimental I-V curves. Figure 1a shows the result of such fitting. Structures for the comparison were MBE grown and processed into real devices, and for this study the samples were selected which give overlapped I-V curves for different waveguide cross-sections what ensures that leakage currents are minimized. The values D = 0.19 nm and K = 9 nm used for best fitting are quite reliable and close to independent estimates (Franckié et al. 2015). Reasonable choice of D and K is further confirmed by the calculated gain spectrum gðhmÞ (see Fig. 1b) which (i) peaks at hm ¼ 0:127 meV, corresponding very well to the experimental wavelength of k ¼ 9:6 lm, and (ii) has the FWHM of 13 meV in close agreement with the measured FWHM = 12 meV of electro-luminescence (Page et al. 2001). The value of the gain peak can be traced as a function of bias. As shown in Fig. 1c, it reaches the value of overall losses divided by the confinement factor ða m þ a w Þ=C ffi 60 cm À1 (Sirtori et al. 1999) for J ¼ J th ¼ 3:61 kA=cm 2 which gives good estimate of experimental threshold current at 77 K (see Fig. 1c). Figure 1d demonstrates further that the agreement between calculated and experimental threshold currents is very good in a wide range of temperatures.
Experimentally inaccessible quantities obtained in the simulations are shown in Figs. 2 and 3. With such data, the in-depth analysis of the electronic transport can be made. In the position-resolved DOE calculated for U = 192 mV/period [overall voltage is 36 (cascade no.) 9 U = 6.91 V], shown in Fig. 2a, the occupation n 3 of upper state 3 is visible as opposed to the occupation of the lower state 2. This means that population inversion (and gain) exists for this bias. Figure 2b reveals that while the carriers in the upper state (3) are thermalized, the distribution of the carriers in the lower state (2) is highly nonthermal: High-k states are largely populated and phonon spectrum (local maxima spaced by DE ¼ hx LO ) is clearly visible. It can be hardly described by the Fermi-Dirac statistics-an approach commonly used in the simplified descriptions-even if the concept of the electronic temperature higher than the lattice temperature is assumed for the laser subbands. The distributions of the currents crossing injection or extraction barriers over E À k 2 space are shown in Fig. 3a and c. From Fig. 3c one can learn that carriers can leave active region in non-thermalized state while from Fig. 3a that they populate upper state being well-termalized. Then, thermalization must take place in the injector wells. This is clearly visible in Fig. 3b which shows position-resolved energetic spectrum of current flowing through the structure. This figure demonstrates also that there is no leakage through excited state no. 4 (which was postulated), and so the injection efficiency to state 3 is close to unity, g 3 ffi 1.
Data in Figs. 2 and 3 enable the estimation of upper state lifetime s 3 which is an important QCL parameter. s 3 can be found from the relation J ¼ en 3 =s 3 plotting upper state population n 3 versus current density J (Fig. 2c). The lifetime s 3 ¼ 0:66 ps we get is significantly shorter than the value s 3 ¼ 1:4 ps calculated from the form factors and assumed in the designing process (Page et al. 2001). In spite of this, as shown in Fig. 2a findings, this result demonstrates that part of the assumptions made in rate-equation analysis is not necessary for the QCLs to work.