Light-enhanced incoherence of electronic transport in quantum cascade lasers

Since their invention in the middle of the 1990s, quantum cascade lasers (QCLs) attract increasing theoretical interest stimulated by their widening applications. One of the key theoretical issues is the optimization of electronic transport which in most of these devices is governed by the injection barrier of QCL heterostructure. In the paper, the nonequilibrium Green’s function formalism is used to study electronic transition through the injection barrier as a function of laser field in the cavity; for the increasing field, a crossover is observed from the strong coupling regime, in which electronic transport through the barrier is coherent, to the weak coupling regime, in which electronic transport gets incoherent. This crossover is characterized by gain recovery time, τrec, which takes sub-picosecond values for mid-IR QCLs operating at room temperature. This time is also important for the performance of devices under steady-state conditions; the maximum output power is obtained when the figure of merit, FOM = (g(0)/gth − 1)/gcτrec [g(0) is the linear response gain, gth is the threshold gain needed to compensate all losses, gc is the gain cross-section], reaches maximum. It is shown that the use of this optimization criterion can result in the structures essentially different from those which can be obtained when the optimized quantity is the linear response gain, g(0).

Quantum cascade laser (QCL) is a semiconductor device used to emit the radiation in the mid-IR or THz frequency range. Opposed to interband devices, it relays on the intersubband transitions occurring exclusively in the conduction band 1 . The population inversion between the laser subbands, which is necessary for the laser action, arises due to the accumulation of electrons in the upper laser subband and very fast extraction of these electrons from the lower laser subband. The upper subband is fed by the electrons tunneling resonantly from the injector through the injection barrier. The QCL emitting in the mid-IR range usually works in the so-called strong-coupling regime 2 . This description refers to the coupling of electronic states located on the opposite sides of the injection barrier; namely: the injector ground state and the upper laser state. When these states are strongly coupled, the carrier transport through the injection barrier is coherent; the nonstationary superposition of the wavefunctions of the states involved in the tunneling transition forms a wave-packet which performs several oscillations before its phase is lost due to the inelastic out-scattering from the upper laser state [3][4][5][6] . Then, the transition through the barrier can be considered as low-impedance if compared to the upper-to-lower transition caused by the inelastic out-scattering. Consequently, both states involved in the tunneling transition are approximately equally populated, and the current is limited by the inelastic scattering 2 . This is equivalent to the assumptions made in the semi-classical description of electronic transport 7,8 . Therefore, it is believed that the methods relying on the Boltzmann transport equation, like rate equation (RE) or Monte Carlo (MC), are sufficient for realistic modeling of mid-IR QCLs 3,7 . The situation may change when the interaction with electromagnetic (e-m) field occurs. For the wavelengths corresponding to lasing transition, the radiative out-scatterings from the upper state become more frequent and may push the device into the weak-coupling regime where electronic transport is incoherent. In this regime, electrons tunneling from the injector lose their phase before they reach the upper lasing state, and so the tunneling transition becomes a limiting factor. The maximum current that can flow in the structure is reduced, and so the optical power. To prevent this situation, the injection barrier is optimized in such a way that the device remains in the strong coupling regime, even if it lases. The optimization of the barrier is not an easy task because it implies the use of the model that works well both in weak and strong coupling regimes. Examples that use the density matrix (DM) method were proposed by Dupont et al. 9 for THz QCL and by Khurgin et al. 10 and Dinh et al. 11 for mid-IR devices. In these attempts, the optimization of the injection barrier aims at the maximization of the gain peak value. This value can be reasonably estimated only when a method that can account for the quantum aspects of electronic transport is used. DM is one of such methods; however, it needs some phenomenological dephasing times which -to some extent -are arbitrarily assumed 12 . Then, the results of the optimization made with DM method suffer from this arbitrariness.
In this paper, another method, that equally treats tunneling and scattering, is employed to study the crossover between the coupling regimes in mid-IR QCL forced by the interaction with the e-m field. This method uses the nonequilibrium Green's functions (NEGF) theory which was formulated in the 1960s 13,14 , and since then has been successfully applied to many transport problems, including the ones in QCLs [15][16][17] . When using the NEGF method, the phase-breaking time does not need to be arbitrarily assumed, or speculated from other considerations, because it inherently arises as a result of the interplay between coherent (quantum tunneling) and incoherent (scattering) processes. Quantities related to this time, like linewidth of gain/absorption spectrum, are then reasonably estimated and do not depend on arbitrarily assumed parameters. In this study, the NEGF method is applied to the number of real QCL structures emitting light of ~5.2 μm wavelength that were tailored for the use in the systems of nitric oxide detection. It is shown that some of these, otherwise very efficient devices, suffer from potential incoherence latent in their structures.

QCL modeling
Density matrix model. The issues raised in the Introduction can be illustrated with a simple 3-state density matrix model, in which the injector state (no. 1) with energy E 1 and upper laser state (no. 3) with energy E 3 are resonantly coupled by the injection barrier. The strength of this coupling is given by the coupling energy Ω 13 . The carriers are out-scattered from the upper state either to the lower laser state (no. 2) with the rate 1/τ 32 or directly to the injector state of the next stage of the cascade with the rate 1/τ 31 . State 1 can also be populated from state 2 with the rate 1/τ 21 . Apart from these non-radiative channels, a radiative channel for → 3 2 transition can be opened by the laser field in the cavity. If the population inversion, Δ = − n n n 32 3 2 , between the populations n i in states 3 and 2 goes to positive values, the flux of electrons flowing in this channel is given by the product of gain cross-section g c , photon flux density Φ (per unit active region width), length d of one QCL module, and Δn 32 . The first three define the field-dependent stimulated emission rate, τ = Φ g d 1/ c st . In this model, the spontaneous emission (much lower than the stimulated one) is not accounted for. When periodic boundary conditions are imposed and the total population in all states is normalized to unity, the steady-state particle current flowing through such system is given by where ρ ij are the elements of the density matrix, τ tr is the effective transport time, and is the tunneling time 18 . In Eq. (2), τ || denotes the dephasing time at the injection barrier, and Δ = − E E ( ) / 1 3  is the detuning from the resonance of 1-3 doublet. The mentioned arbitrariness of the DM model is stored in the time τ || ; τ −1 is half the sum of the state decaying rates due to intersubband processes plus some "pure" dephasing rate, which is a phenomenological component 4 . Substituting Eq. (2) into Eq. (1), one gets the well-known Kazarinov-Suris relation 19 in a form adapted for the current 3-state system For the lasing device, the amount of the current increment, Φ − Φ = j j ( ) ( 0), caused by the laser field is of interest, because it is directly related to the flux of particles that can emit the photon and so to the output power. In Fig. 1, this quantity is plotted using Eqs. , reasonable for mid-IR QCL. The calculations were done at the resonance, i.e., for Δ = 0, so the plots represent the maximum increase of the current and the maximum optical power available in the system.
The coupling regimes are defined by the relation between the terms in the denominator of Eq. (3) 2 . Namely, the strong coupling (coherent transport) occurs when the last term largely exceeds the first two. Then, initially (i.e., for , the system is in the coherent regime for both values of the coupling energy, although for ħΩ 13 = 2 meV the system is "less coherent", i.e., it has lower value of  τ τ Ω * 13 2 3 . In the strong coupling regime: the current j is exclusively described by the scattering lifetime τ ⁎ 3 , the occupations ρ = n ( ) i ii of states 1 and 3 are approximately equal and -as illustrated in Fig. 2a -electrons perform several oscillations before they settle at the steady-state value. As shown in Fig. 2b, the coherence is destroyed by the laser field; in this case, the oscillations are more damped and disappear after one cycle. As demonstrated in Fig. 1, the current increases with the laser field; however, the increase of the current is lower for initially less coherent system than for the other one. Worth to note is that both oscillatory and non-oscillatory responses were observed in the number of pump-probe experiments 5 , also depends on this time, which one would like to avoid as its value comprises phenomenological components. The NEGF model presented in the next section does not have such disadvantage. NEGF model. In the NEGF approach, the dephasing time does not appear per se. States and gain broadenings arise exclusively from all inelastic and elastic scatterings included in the model as well as due to the fact that the quantum system is open. All these factors are taken into account through theoretically well-justified physical models: the system Hamiltonian and the scattering selfenergies. As the aim is to derive conclusions for real QCL structures, the NEGF model of the device presented in this section is much more realistic than the simple 3-state system considered in the previous section.
So far, all QCL nanostructures are n-type. They consist of parallel layers made of different materials. The appropriate Hamiltonian is then one-band 1D operator in the growth direction z, parametrized for the in-plane momentum modulus, = k k xy 1,16 . In Eq. (4), the potential V z ( ) includes both the variation of the conduction band edge E c (z) and the Hartree potential. Nonparabolicity is accounted for through the energy-dependent effective mass,  is the bandgap energy. Other symbols have the usual meaning. It was shown that such Hamiltonian gives results that very well match the results obtained with the eight-band kp model 1 .
The equations of the NEGF formalism 13,14 applied to the Hamiltonian of Eq. (4) were solved numerically in the real space (position basis) for the steady-state. The scattering selfenergies in the formulations, suitable for the Hamiltonian of Eq. (4), were provided in the number of papers 16,22,23 . They cover most important scattering processes involving: optical and acoustic phonons, ionized impurities, alloy disorder, and interface roughness. Recently, the selfenergies for the interaction with the optical field have been also provided 24 , what enables this study. As already mentioned, the electron-electron interactions were treated within the mean-field approximation. This was done by solving the Poisson equation self-consistently with the NEGF equations. For the former, the boundary conditions ensured device neutrality. At the output of the procedure, the 4-argument Green's functions ′ G z z E k ( , , , ) were obtained. Physical quantities were calculated from functions G R and G < . Appropriate formulas for the densities of currents, states, and carriers are provided, e.g., in ref. 16 . The linear response (small signal) gain/absorption was calculated based on the theory outlined in ref. 15 , adapted for the case of energy-dependent effective mass 25 . More details of NEGF calculations are provided in the Methods section.
Devices. The attention is paid to two QCL structures that emit IR radiation at ~5.2 μm wavelength. Namely: structure D of Diehl et al. 26 , and structure E of Evans et al. 27 . Both structures are very efficient. However, the output power of ~0.3 W emitted at room temperature from structure D is only approximately half of that emitted from structure E under similar operating conditions. We claim that the poorer performance of structure D is due to the larger incoherence of tunneling transition forced by the laser field in this structure. To prove this claim, the simulations were performed; however not for original devices. For those used in the simulations, all parameters, which can influence the output power, were kept identical. Consequently, the calculations were performed for devices which have equal number of periods, , was assumed for both devices. This value corresponds to overall losses: α m = 2 cm −1 , α w = 3.8 cm −1 in a typical waveguide, with a confinement factor, Γ = . 0 6 26 . Under such assumptions, the output optical power was calculated as 28 p th where w = 7.5 μm is the cavity width, R = 0.27 is the facet reflectivity, and E γ is the photon energy. As d and E γ are very similar for both structures, the output power depends exclusively on the value of the photon flux Φ th required to clamp the gain to its threshold value g th . In general, Φ th can be different for different structures.

Results and Discussion
Results of the simulations presented in Fig. 3a confirm the experimental observation; the output power available from structure D is lower by 30-40% than from structure E. The absolute values of the optical power are also not very far from those reported in the experiment, although the differences are expected; real devices have different waveguide dimensions and are doped to different densities. They work continuous wave (cw), and so suffer from self-heating. On the contrary, the simulations assume that the dissipated electrical power does not heat up the structures. www.nature.com/scientificreports www.nature.com/scientificreports/ Data in Fig. 3 illustrate how the light-current-voltage characteristics were obtained. Initially, current-voltage (I-V) and gain-current characteristics were calculated for the laser field turned off. In the second step, the laser field was turned on and increased until the gain peak was clamped to the threshold value, = .
− g 9 5 cm th 1 . Data for the "off " state are also included in Fig. 3. One can observe that the amount of excess current, I on -I off , that can be converted into photon flux is larger for structure E than for structure D. According to the DM model, this can be attributed to more coherent transport in structure E. Indeed, the calculations of the coupling energy performed like in ref. 29 give the values: Ω = .
3 55 meV 13  for structure E. The difference is due to different thicknesses of the injection barrier, which is 4.2 nm for structure D and 2.9 nm for structure E.
One can observe in Fig. 3b that for both structures the gain peak in the "off " state reaches approximately the same value of 30 cm −1 . While lasing, the gain is clamped to the value g th , which is also equal for both structures. As the increase of the current, I on -I off , is lower for structure D, the suppression of the gain by the optical field must be more "efficient" in this structure. This "efficiency" can be quantified by the gain recovery time, τ rec , appearing in the relation 30 c rec or alternatively by the saturation intensity, Φ s = 1/(τ rec g c d) 31 . According to the DM model  30 . For the structures in question, the gain recovery time can be evaluated from the data presented in Fig. 4. In the simulations made for this purpose, the bias voltage was kept at a constant value (in order to fix the tunneling time), while the photon flux was increased and gain spectrum was calculated for certain values of the photon flux. In general, the dependence observed in Fig. 4 demonstrates the decrease of the gain caused by the optical field. As can be seen, the gain suppression is stronger for structure D, as expected for "less coherent" structure. The curves differ and larger flux is needed in order to clamp the gain to the threshold value for structure E. Therefore, this structure emits more power. Gain recovery time estimated by fitting the data in Fig. 4 to Eq. (6) takes the value τ = .
0 34 ps rec for structure D and only τ = . 0 24 ps rec for structure E. The values of gain cross-section necessary to get these estimates were obtained in an independent numerical simulation described in the Methods section. The estimates of τ rec are a bit lower than the value 0.47 ps found in a similar numerical experiment for another mid-IR QCL 32 . The larger value of the gain recovery time implies that the structure considered in ref. 32 was even "more incoherent". Indeed, the value τ τ Ω ≅ . * 0 2 13 2 3  estimated for this structure 32 classifies it deeply in the weak coupling regime. The observed dependency: the stronger coupling, the shorter gain recovery time suggests that the latter cannot be exclusively connected with the scattering times as predicted by semi-classical models 30,32 . As stems directly from Eq. (7), this time also depends on the tunneling time which is related to the strength of the coupling energy, Ω 13 , through Eq. (2). www.nature.com/scientificreports www.nature.com/scientificreports/

Structure optimization
The analysis performed in the previous section proves that the difference in the performance of structures D and E results from the degree of coherence (of tunneling injector-to-upper-state transition) latent in these structures. The measure of this performance is the gain recovery time which should be minimized to obtain more optical power. The correct strategy for structures optimization (if the output power is a criterion) is not only to maximize the gain (like in refs. [9][10][11] ) but simultaneously to minimize τ rec . Quantitatively, the best method for this would be to maximize the quantity , which corresponds to the very thin/ low barrier. In this case, the gain spectrum broadens at the expense of the gain peak 33 . This effect is not accounted for in the DM model. To take it into account, a more detailed quantum model of electronic transport, either DM 29 or NEGF, is necessary. A DM approach that used gain peak as a criterion (but not Eq. (8)) was presented, e.g., in refs. [9][10][11] . In the NEGF method, the value of Φ th can be directly calculated (see Fig. 4) without the need of assuming any phenomenological parameters or even the estimation of the quantities, like τ rec or g c . The results of calculations that maximize the figure of merit of Eq. (8) for structure D with different barrier thicknesses are presented in Fig. 5. As can be seen, Φ th goes through a broad maximum which means that finding optimal barrier thickness is possible. Indeed, the maximum output power P = 0.62 W found for structures with 3.2 nm-thick barrier is larger than the maximum output power found for structure with 4.2 nm (0.46 W) and 2.4 nm (0.52 W) thick barriers (see Fig. 6). In Fig. 5, the results of optimization, when the linear response gain g(0) is the criterion, are also included. As can be seen, they may lead to different conclusions concerning optimal barrier thickness.

conclusions
The dependence shown in Fig. 4 demonstrates that the implementation of light-matter interaction in NEGF formalism through appropriate selfenergies enables going beyond linear response to an external e-m field without the need of solving the time-dependent Kadanoff-Baym equation. While the "selfenergy" approach has already been used for interband devices [34][35][36] , it is new for lasers relying on intersubband transitions. Till now, the light-matter interaction was implemented in the NEGF method applied to QCLs through the decomposition of the Green's functions and selfenergies into higher harmonics of the lasing frequency 30,32,37 . This approach is more exact, but generates much larger numerical load than the approach presented in the current study.
It was shown that the laser action pushes the structure into the weak coupling regime where transport through the injection barrier is less coherent and limited by tunneling. Numerical simulations provide unique tool to study this light-induced crossover; the upper state lifetime which is a combination of non-radiative and radiative scattering rates can be tuned by changing the photon flux, keeping the remaining parameters unchanged. In working devices, such tuning is hardly possible, because photon flux saturates at the value which clamps the gain to its threshold value. The simulations show that light-enhanced incoherence of the tunneling transition is characterized by the gain recovery time which increases with the decreasing initial (static) coherence latent in the quantum design of the structure. Contrary to semi-classical results, this time depends not only on the scattering lifetimes but also on the tunneling time through the injection barrier. Due to this dependence, the devices (even those emitting in the mid-IR range) may experience the limitations introduced by too long gain recovery time. Real devices suffering from this limitation do exist 26 . To increase their performance, the strategy for barrier optimization should be reconsidered; the criterion of maximizing the linear response gain (gain in the no-lasing state) should be replaced with maximizing the figure of merit provided by Eq. (8). Data presented in Fig. 5 illustrate that − g 9 5 cm th 1 (circles -left axis) and gain peak at zero flux g(0) (triangles -right axis), calculated with the NEGF method for structure D with various injection barrier thicknesses. Structures were biased with the voltage U = 320 mV/period. these two criteria can lead to different optimal structures. Arriving at these conclusions was possible due to full quantum studies of QCL nanostructures which were done with NEGF method.

Methods
Results for the DM model were obtained by the analytical solution of the Liouville equation completed with the terms that characterize scattering and dephasing of density matrix elements For the 3-state model, discussed in the paper, H, ρ and S are 3 × 3 matrices.
The equations of the NEGF formalism in the position basis were adapted from ref. 16 . They were applied to the structures truncated to a bit more than one QCL module (period). The rest of the cascade was mimicked by suitable contact selfenergies 17 . Non-uniform discretization mesh was used to reduce numerical load without losing the accuracy 38 . The parameters m* and E g for ternary materials used for QCL layers were taken from ref. 39 . Conduction bandgap offsets were calculated using model solid theory including strain 40 and material parameters that were taken from ref. 41 . This method was also used to evaluate conduction band offsets used in the parametrization of alloy disorder scattering. The scattering selfenergies were implemented in a way described in refs. 22,23 for the scatterings from ionized impurities, alloy disorder, LO-phonon, and interface roughness. For scattering from  www.nature.com/scientificreports www.nature.com/scientificreports/ acoustic phonons, the energy-averaged approximation of ref. 16 was implemented. The electron-photon selfenergies were calculated using the low-density approximation 24 . All parameters used in the NEGF calculations, as well as the set of microscopic results obtained with the NEGF method for structure D, are provided in the Supplementary Information.
The gain cross-section at the lasing energy E γ can be estimated as the ratio of the gain at this energy and the population inversion, = ∆ γ γ g E g E n ( ) ( )/ c 32 . This definition may be applied to the total subband populations and total gain only when bands are parabolic and the populations are thermalized. In mid-IR QCLs, energy levels are well above the conduction band edge and nonparabolicity must be taken into account. It was also shown that the populations in QCL subbands are highly nonthermal 42 . In effect, not all carriers contributing to Δn 32 contribute to g(E γ ). To get rid of such particles, g c was estimated from momentum-resolved calculations, = ∆ γ γ g E k g E k n k ( , ) ( , )/ ( ) c 32 . The results of these calculations are presented in Fig. 7. Except for the largest momenta, g c takes an almost momentum-independent value of ~0.60 × 10 8 cm, quite close to the value 0.57 × 10 −8 cm found for similar (GaInAs/AlInAs) QCL 28 . At largest k's, g c drops to lower and negative values since in this range: (i) the gain at the lasing frequency E γ is not produced (because = − < γ E E E E 32 3 2 due to bands nonparabolicity); (ii) the populations are not inverted (see Fig. 7a). For the comparison; if the total (i.e., summed over all k) populations in the laser subbands and total gain were used, the (largely) overestimated value, g c = 1.53 × 10 −8 cm, was obtained.