Coupled carrier–phonon nonequilibrium dynamics in terahertz quantum cascade lasers: a Monte Carlo analysis

The operation of state-of-the-art optoelectronic quantum devices may be significantly affected by the presence of a nonequilibrium quasiparticle population to which the carrier subsystem is unavoidably coupled. This situation is particularly evident in new-generation semiconductor-heterostructure-based quantum emitters, operating both in the mid-infrared as well as in the terahertz (THz) region of the electromagnetic spectrum. In this paper, we present a Monte Carlo-based global kinetic approach, suitable for the investigation of a combined carrier–phonon nonequilibrium dynamics in realistic devices, and discuss its application with a prototypical resonant-phonon THz emitting quantum cascade laser design.


Introduction
Within the conventional treatment of carrier-quasiparticle interaction in optoelectronic devices, the typical assumption is to consider the quasiparticle subsystem (e.g. photons, phonons, etc) as characterized by a huge number of degrees of freedom if compared to those of the carrier subsystem. In other words, the former behaves as a thermal bath: it is always in thermal equilibrium and is not significantly perturbed by the carrier dynamics. For many ultrafast and/or high-field-transport phenomena in nanomaterials and nanodevices, however, this assumption is highly questionable, and the carrier-quasiparticle treatment should be extended to account for quasiparticle nonequilibrium regimes as well [1].
Generally speaking, the latter may be classified according to the following qualitative subdivision. On one hand, we deal with coherent quasiparticle phenomena, where the quasiparticle non-diagonal density-matrix elements play a crucial role; in addition to various coherent-light effects, as well as to quantum-optics phenomena, a typical example is that of the so-called coherent phonons. On the other hand, also in the absence of such phase coherence effects, the quasiparticle diagonal density-matrix terms may deviate significantly from the equilibrium Bose distribution; two relevant examples are photonic cavity-mode amplification and hot-phonon effects.
Quantum cascade lasers (QCLs) offer a paradigmatic behaviour within the second class of the above-mentioned nonequilibrium phenomena. In the core region of these devices, the dynamics of the charge carriers-and therefore the light emission-is not only unavoidably affected but also intentionally controlled by exploiting the electron-phonon interaction. This is particularly evident for mid-infrared emitting designs [2], but is a possible strategy also for terahertz (THz) operating devices, where the proximity of the photon and longitudinal-optical (LO) phonon energy initially motivated the search for alternative, and, in some sense, less intuitive, schemes [3,4].
The principle of operation of QCLs exploits downhill electronic transitions along the biased, suitably designed ladder of subbands, in which the core-heterostructure conduction band is split due to quantum confinement. In particular, the gain regime in the active region is established and maintained by means of a properly tailored electron energy relaxation dynamics, resulting in a selective depopulation of the diverse subband states. In conventional mid-infrared as well as in THz resonant-phonon designs, the separation between the lower laser and the ground subband within each period matches the LO-phonon energy so as to maximize carrier relaxation out of the former into the latter via LO-phonon emission. In this case, a significant nonequilibrium population of LO phonons is present in the core region of the operating device when, as is generally the case, their generation rate is larger than the anharmonic decay rate into acoustic modes. Indeed, by comparing the potential energy drop required for the proper subband alignment with the LO-phonon energy, it appears that, while cascading along the sequence of periodically repeated identical stages within the core region of resonant-phonon THz QCL designs [5], each electron generates at least one optical phonon per period. Since the state-of-the-art designs may include from several tens up to more than one hundred stages, it appears that these devices, besides being reliable and robust infrared light sources, are also rather efficient LO-phonon generators. From the theoretical point of view, the combined effect of high LO-phonon generation rates and limited thermal conductivity makes it imperative to develop a realistic description of charge transport in these devices that takes into account the carrier interaction with a nonequilibrium phonon subsystem.
The impact of nonequilibrium phonon 2 effects on the electron relaxation dynamics, and therefore on the performances of QCLs, has been theoretically addressed in the past, in terms of Monte Carlo (MC) kinetic treatments limited to a subset of subbands [6][7][8], within a prototypical Krönig-Penney model [7,8] or via partially refined implementations [9]. Recent experimental studies have pointed out distinctive nonequilibrium phonon features both in mid-infrared [10] as well as in THz QCLs [11,12], evidencing a significant heating of both the electronic and lattice systems [13]. Motivated by such experimental findings, in this paper, we will present a MC-based global kinetic theoretical approach that we have developed for the analysis of these effects in realistic devices [14,15]; the latter includes on the same footing the fully threedimensional nature of the multisubband electronic structure as well as of the phonon degrees of freedom, and is able to overcome significant limitations affecting some of the MC techniques previously mentioned.
The paper is organized as follows. In section 2, we shall describe the theoretical framework for the microscopic modelling of the coupled carrier-quasiparticle dynamics in semiconductor nanostructures, based on a closed set of equations for the corresponding distributions; in section 3, we shall apply our approach to a prototypical THz emitter; finally, in section 4, we shall summarize and draw a few concluding remarks.

Physical system and theoretical description
The description of the coupled-phonon nonequilibrium dynamics in our prototypical THz emitting QCL may be performed within a purely semiclassical picture for both subsystems, in which the relevant kinetic variables are the electron and phonon distributions. In particular, a coupled set of nonlinear equations of motion for the latter may be derived following the conventional single-particle correlation expansion of the quantum-kinetic theory, and neglecting all non-diagonal carrier as well as quasiparticle density-matrix elements [16]. Actually, while coherent phonons and carrier phase coherence effects may be investigated at the transient ultrafast timescale [17], they are expected to play a minor role in the steady-state regime we are presently interested in. Indeed, in spite of the fact that the system under investigation is spatially inhomogeneous, the LO-phonon subsystem can still be described via a diagonal density matrix in q, parameterized by a meso/macroscopic spatial coordinate related to the thermal-transport space scale. This scenario is similar to the conventional semiclassical or Boltzmann transport picture, where the carrier distribution function and its scattering dynamics can be parameterized by a real-space coordinate, but such nonhomogeneities are characterized by a space scale much longer than the carrier coherence length.
In principle, a proper treatment of the steady-state properties of our device should then start from both the electron and the phonon coupled Boltzmann equations. This task is quite demanding, since it requires, for the latter, the inclusion of both acoustic and optical modes, preferably with finite size and quantization effects. To start focusing on the main physical aspects of the energy redistribution between charge carriers and lattice degrees of freedom, in this paper we consider the full electron subsystem (i.e. the complete set of active region and injector subbands within each period) interacting with bulk LO-phonon modes, whose decay into acoustic modes (lattice thermal bath at temperature T L ) is accounted for by means of a phenomenological lifetime τ .
In other words, energy dissipation occurs-and is modelled-as a two-step process: the carrier subsystem transfers a significant amount of energy to the quasiparticle (LO phonons) one; the latter, in turn, will transfer its excess energy to additional degrees of freedom (acoustic phonons) that are not taken into account dynamically and which are characterized by a much larger heat capacity. The latter phenomena are described via a standard relaxation-time model. The structure of the closed set of coupled equations corresponding to the above-described dynamics is then the following: where f νk is the electron distribution function corresponding to the single-particle state in subband ν and with in-plane wavevector k , while n q is the average phonon occupation number pertinent to a single LO-phonon mode with wavevector q.
Within the Fermi's golden rule approximation, the carrier-phonon coupling Hamiltonian (e-LO) translates into the typical Boltzmann scattering term, consisting of in-and out-scattering contributions, due to phonon emission and absorption processes. Focusing first on (1), the latter have the form Here where g ± νk ,ν k ;q are the matrix elements (within our multiband electronic basis νk ) of the (q-dependent) carrier-LO phonon coupling coefficients, and the sign + (−) refers to emission (absorption) processes.
The remaining term, labelled with s, at the right-hand side of (1), accounts for further scattering mechanisms, of both intrinsic and extrinsic character, that may affect the carrier dynamics. In particular, within the former type of processes, the electron-electron interaction has proven to generally play an important role in these devices [3] and will therefore be included in our present analysis. Conversely, extrinsic mechanisms, such as interface roughness and carrier-impurity scattering, should be considered if one is interested in a strictly quantitative analysis of a specific device behaviour. Such processes, indeed, do not significantly modify the trend of the current-voltage characteristics since, opposite to carrier-LO phonon scattering, they have a threshold-less nature and are poorly dependent on the applied bias. Moreover, they are strongly device/sample dependent and would therefore require a phenomenological description. For these reasons, they will not be considered in the remainder of this paper.
The structure of equation (1) allows for a self-consistent charge-conserving MC carrier transport simulation, provided that the phonon distribution n q is known [17]. A typical example is when the latter may be reasonably well approximated by the constant-i.e. q independent-value of the Bose-Einstein distribution at a given (quasi)equilibrium temperature. This is, however, not the case when the interplay between carrier and phonon dynamics is such that the quasiparticle distribution is significantly driven out of equilibrium during device operation. If nonequilibrium phonons are considered, their occupation numbers n q are no longer q independent and have to be evaluated by solving the corresponding dynamical equation (2).
The Boltzmann transport equation (1) for the electron subsystem, with the contribution in (3), goes then together with the phonon counterpart in (2). The explicit form of the key quantities entering the latter is Here, the − n q −n th q τ term accounts for loss processes from the LO-phonon subsystem due to anharmonic decay into acoustic modes and n th q being the Bose-Einstein distribution at the lattice (acoustic) bath temperature T L .
The structure of the above-described coupled-equation formalism is suitable for a MC particle-like approach to both the electron and phonon dynamical equations. Indeed, the latter may be conveniently sampled by means of a generalized MC simulation able to deal with a fixed number of electrons and a variable number of phonons, and therefore implemented in the stateof-the-art simulation tools [17,18]. In particular, the q-dependence of the phonon occupation numbers can be accounted for by means of a combined self-scattering and rejection technique. Moreover, the use of a periodic boundary condition approach allows us to evaluate the device performances in a 'closed-circuit' scheme, in which the only free parameters are the LO-phonon lifetime τ and the lattice temperature T L . The latter may be accessed by means of state-of-theart microprobe band-to-band photoluminescence experiments [12,13], while realistic values for the former in bulk materials are in the range 6-9 ps [19].

Application to a prototypical terahertz emitter
During the operation of a typical QCL structure, a significant population of LO phonons is generated and a considerable amount of energy is dissipated in the device Joule heating. This is, in particular, evident in resonant-phonon designs, where the threshold nature of the LOphonon emission process is intentionally exploited to selectively depopulate the lower laser subband. In these devices, charge transport is mainly governed by carrier-phonon scattering and the gain can be principally attributed to a single transition within each period. Indeed, the combination of high LO-phonon generation rates and limited thermal conductivity produces a relevant nonequilibrium LO-phonon distribution that unavoidably affects the electron relaxation kinetics and therefore the device performances.
To have a deeper and more quantitative insight into the impact of 'hot' phonon phenomena on the nonequilibrium electron dynamics in THz QCL emitters, we have applied the microscopic modelling framework presented in section 2 to the resonant phonon device based on the quantum design described in [14,15,20]. As anticipated, besides carrier-LO phonon scattering, we have included also carrier-carrier interaction in the electron dynamical equation (1), employing the well-established time-dependent static-screening model commonly adopted in two-dimensional systems [21]. The lifetime τ in (2) is set to the value of 6 ps, constant throughout the whole device operation range. As far as the remaining input parameter, T L , is concerned, its values have been taken directly from the experimental characterization of the operating sample by means of microprobe band-to-band photoluminescence; the latter evidences a significant heating of the device, with T L being bias dependent and locally higher than the cryostat temperature (100-180 K-varying with the applied bias-versus 80 K) [14].
An important feature of the present-global and three dimensional-kinetic approach is the ability to directly access the electron distribution functions corresponding to the diverse subbands as a pure output of the MC simulation. In other words, no assumption of intrasubband thermalization is a priori made. The latter, instead, may eventually show up due to the synergistic interplay between carrier-carrier and carrier-LO phonon scattering. This is indeed the case also for the device we are considering, where the simulated distribution functions for the various subbands exhibit a typical heated Maxwellian form, such as the one reported in figure 1. A proper fit of the latter allows one to estimate the effective temperature characterizing the electron energy distribution within each subband. If this analysis is systematically performed, it turns out that a 'hot electron' regime is established in the operating device, in which the carrier effective temperatures are higher than the lattice one; this theoretical result is in good agreement with experimental findings [14].
While the above-described heating of the electronic subsystem may be qualitatively expected, its quantitative impact does strongly depend on the device under examination, being the result of a non-trivial interplay between single-particle (carrier-LO phonon) and two-particle (carrier-carrier) interactions. In particular, the efficiency of the latter mechanism in setting up a heated distribution is strongly density dependent and its effect cannot be a priori quantified. In other words, a 'hot' (in the broader sense) carrier subsystem is a reasonable expectation, but its thermalized nature is not so straightforward and can only be confirmed via a kinetic treatment like the one employed in this paper. In this case, the diverse temperatures characterizing the quasi-equilibrium distributions within the various subbands may also be extracted from the MC simulation: the latter are key ingredients in any simplified n-level (i.e. one-dimensional) reduced description of the three-dimensional device. As will be discussed later in the paper, a quite different scenario appears on the side of the LO-phonon subsystem, which presents a definitely out-of-equilibrium behaviour and cannot be parameterized in terms of any effective temperature.
In addition to the carrier heating, our MC simulations also show that a population inversion regime is found between the upper and the lower laser subbands throughout the considered bias range, again in agreement with experimental evidence. A key question that recurrently occurs when investigating novel QCL designs regards the feedback that the nonequilibrium phonon dynamics produces on the device gain performances. An unavoidable reduction of the latter is indeed expected, but its analysis and quantitative estimate require all the power and flexibility of our MC-based approach, which allows one to selectively switch on and off the diverse contributions to the dynamical equations, highlighting their individual and synergistic roles. Coming to the THz emitting device considered here, indeed we find that the population inversion is reduced to approximately 85% of the value corresponding to the ideal case, in which the electron subsystem interacts with an equilibrium LO-phonon population at the cryostat temperature. A closer look at the subband population dynamics shows that this is the consequence of a reduced extraction efficiency out of the lower laser subband (thermal backfilling) into the ground injector one, rather than of the thermally activated depopulation of the upper laser states. In other words, the presence of a 'hot' phonon population leads to an enhancement of the phonon absorption rates, thereby reducing the net energy relaxation/loss for the charge carriers [22].
The simulated current-voltage characteristics of our prototypical device are shown in figure 2. In the latter, indeed, results obtained at diverse levels of description are compared: one set of data (triangles) corresponds to the case, in which the complete (acoustic plus optical) phonon subsystem behaves as a thermal bath at the cryostat (80 K) temperature, while the other two describe the charge transport response when LO-phonon nonequilibrium phenomena are dynamically taken into account and the acoustic bath is at equilibrium either at the cryostat temperature (squares) or at the higher (measured) ones (discs). All the curves exhibit the typical behaviour of QCL structures, i.e. the current density increases on increasing the bias as long as the injector and upper laser subbands are kept aligned, while a negative differential resistance region shows up at higher fields. However, when comparing the disc data set with the other two, one sees that the above-discussed diminished cooling rate out of the lower laser subband has a beneficial effect from the point of view of the global transport dynamics across the device core region. It in fact prevents the accumulation of charge carriers in the ground injector subband, which is less efficiently coupled to the upper laser one, and forces them to follow the alternative relaxation path into the latter from the other subband of the injector doublet. The resulting current enhancement may amount up to more than a 60% increase with respect to the equilibrium LO-phonon bath case (triangles) and is again a 'hot'-phonon effect in the broader sense. Indeed, the inclusion of nonequilibrium LO phenomena at the bare cryostat temperature for the acoustic phonon bath (squares) only accounts for a small fraction of such correction. Similar to the hot-electron scenario previously mentioned, a closer inspection of the current-voltage profiles reveals once again that the quantitative impact of the hot-phonon dynamics is not only dictated by the specific device design but also depends on the particular operation condition, i.e. on the applied bias. In particular, it follows that it is not possible to employ an effective/heated-phonon-temperature model, since the latter will be strongly device and bias dependent, and thus could only operate as an a posteriori fitting parameter without any predictive added value.
The current density across the device is a global quantity, i.e. it somehow averages over the details of the phonon distribution and does not immediately allow one to discriminate between the effects of a quasi-equilibrium population (thermal, though heated) and a nonthermal one. Indeed, the assessment of a nonequilibrium LO-phonon population can be unambiguously recognized by looking at the details of the solution of equation (2). Figure 3 shows the n q values extracted from the MC simulation, as a function of the in-plane modulus q = |q in plane |, 3 and compares them with the thermal, q independent, Bose distribution at the temperature T L . Although the strictly quantitative values for the former and the latter depend on the applied bias, the general behaviour is common and shows that a significant amount of small-q phonons is generated-and maintained-in our operating device. This is indeed a key feature of the polar phonon emission process, whose rate decays as 1/q 2 therefore privileging, at resonance, q wavevectors close to the zone centre. Such a nonthermal phonon population scenario has been recently probed and confirmed by microRaman spectroscopy analysis, evidencing an extremely good agreement between simulated and measured data [15]. A final remark is mandatory. From a strictly quantitative point of view, the analysis presented so far applies to a specific region of the device, the one that is selectively accessed by means of microprobe experimental techniques. Indeed, the cross-plane thermal conductivity of these multi-interface structures is strongly reduced with respect to the bulk materials and this results in a non-uniform heating along the growth direction. In other words, the measured T L values are local quantities, which may significantly vary in the direction of the applied bias. Since T L is an input parameter, a possible way to account for these thermal resistance effects in the proposed model consists in assuming its value to be position dependent and setting it from an experimental device temperature mapping.

Summary and conclusions
In the present-day optoelectronic quantum devices, excited carriers are often found to emit a large fraction of phonons, thus driving the corresponding distribution out of equilibrium. This is particularly true when relaxation via LO-phonon emission is the mechanism of choice to achieve the desired transport/gain device performances. In this paper, we have analysed hot electron-hot phonon effects in a prototypical THz QCL structure. In particular, we have developed an MC-based global kinetic approach describing the complete interacting electronic subsystem (i.e. the complete set of active-region and injector subbands), coupled to a bulk LO-phonon subsystem, whose decay into acoustic modes is modelled by a phenomenological lifetime. Energy dissipation occurs then in the form of a two-step process: the carrier subsystems first deliver a significant amount of energy into the LO-phonon one; the latter then transfers its excess energy to additional degrees of freedom, characterized by a much larger heat capacity. Simulated results show a very good agreement with experimental data, evidencing how the nonequilibrium phonon population may affect the electro-optical device performances. The proposed scheme may be extended to account for a more refined modelling of the LO-phonon subsystem, and to include a more realistic description of their decay into acoustic phonon modes and of the heat transport dynamics.