Monte Carlo Kinetic Modeling of the Combined Carrier-Phonon Nonequilibrium Dynamics in Semiconductor Heterostructure Devices

Abstract

Electron-phonon interaction is a key mechanism for charge and heat transport in both bulk materials as well as in state-of-the-art electronic and optoelectronic solid-state devices. Indeed, that of an effective heat dissipation, at the diverse design levels, has always been a primary issue in device operation and performances. In various circumstances, the charge carrier subsystem happens to be coupled to a significant nonequilibrium optical phonon population. This regime may be particularly pronounced in new-generation quantum emitters based on semiconductor heterostructures and operating both in the mid-infrared as well as in the terahertz region of the electromagnetic spectrum. In this chapter, we review a global kinetic approach based on a Monte Carlo simulation technique that we have recently proposed for the modeling of the combined carrier-phonon nonequilibrium dynamics in realistic unipolar multisubband device designs. Results for the case of a prototypical resonant-phonon terahertz emitting quantum cascade laser are shown and discussed.

1. Introduction

Electron-phonon scattering is a main mechanism behind the charge carrier dynamics in solid-state materials and devices. Its simplest treatment assumes to consider the phonon subsystem as characterized by a huge number of degrees of freedom in comparison to the electron one. The former, in other words, behaves as a thermal bath, being always in thermal equilibrium and not crucially perturbed by the dynamics of the latter. This premise is however unconvincing for many new-generation devices, involving ultrafast and/or high-field-transport phenomena in nanomaterials and nanostructures. In these cases, the carrier-quasiparticle treatment has to be extended to properly account for phonon nonequilibrium regimes as well [1].

In general terms, the latter may be subdivided into two classes: one including coherent phenomena, where the quasiparticle nondiagonal density matrix elements play a crucial role—a typical example being that of the so-called coherent phonons—, and one in which, even in the absence of such phase-coherence effects, the quasiparticle diagonal density matrix terms significantly deviate from the equilibrium Bose distribution—a typical case being that of the so-called hot phonon effects—.¹

A paradigmatic example within this second class of nonequilibrium phenomena is represented by quantum cascade lasers (QCLs). In the core region of these devices, the carrier dynamics in general, and the population inversion regime in particular, are not only unavoidably affected but in fact intentionally controlled through the electron-phonon interaction. This strategy was the key to success for the pioneering mid-infrared emitting devices [2], but comes out to have a crucial role also for terahertz operating structures, where the closeness of the photon and longitudinal-optical (LO) phonon energy initially oriented toward alternative, and less intuitive, designs [3, 4].

QCLs are unipolar coherent light sources whose basic principle of operation may be declined in a variety of specific—and often rather complex—designs, guided by the requested/desired output and performances [5, 6]. For the purposes of the present review, it will suffice to focus on their core region. The latter is basically a semiconductor multiquantum-well heterostructure consisting of a periodic repetition of identical stages, each made up of an active region sandwiched between an electron-injecting and an electron-collecting region. When a proper bias is applied, each collecting region behaves as an injection one for the subsequent stage, and an “electron cascade” takes place along the energy-level staircase in which the conduction band is split due to quantum confinement. Since their first demonstration, the potentialities of band-gap engineering and molecular beam epitaxy have been inspiring several successful QCL designs for both the active region and the injector/collector part. In particular, the key issue of establishing and maintaining the gain regime in the active region can be accomplished by means of a proper tailoring of electron energy relaxation dynamics, resulting in a selective depopulation of the diverse subband states. Both conventional mid-infrared as well as so-called THz resonant-phonon designs [7] are devised in such a way that the separation between the lower lasing subband and the ground one within each period matches the LO phonon energy. This allows one to maximize carrier relaxation out of the former into the latter via LO phonon emission.

As a by-product, a LO phonon population significantly out of equilibrium is expected within the core region of the running device when—as it often happens—their generation rate is higher than the one describing their anharmonic decay into acoustic phonon modes. Actually, if one compares the LO phonon energy with the potential energy drop allowing for the essential subband alignment along the growth direction, it comes out that, while cascading along the succession of equal stages, each electron generates at least one optical phonon per period. State-of-the-art designs may include from several tens up to more than 100 stages. It therefore easily appears how these devices are not only reliable and robust infrared light sources but also rather effective LO phonon generators; this poses severe heat-removal issues especially in high-performance device architectures [8, 9]. Already at the drawing board stage, it is evident that the combined effect of high LO phonon generation rates and limited thermal conductivity requires a realistic description of charge transport in these devices which goes beyond the thermal bath picture, taking into account the carrier interaction with a nonequilibrium phonon subsystem. In particular, while simplified multilevel rate-equation models may offer a computationally light picture of the electronic cascade in the active region [10], their quantitative/predictive value is limited and does for sure greatly benefit from the inputs offered by the more exhaustive three-dimensional multisubband kinetic picture.

In the past, several theoretical schemes have been proposed to model the impact that nonequilibrium phonon effects may have on the electron relaxation dynamics, and therefore on the performances of QCLs: Monte Carlo (MC) kinetic approaches restricted to a subgroup of subbands [11, 12, 13], within a prototypical Krönig-Penney model [12, 13] or partially refined implementations [14]. From the experimental point of view, various studies have highlighted peculiar nonequilibrium phonon features both in mid-infrared [15] as well as in THz QCLs [16, 17]; in particular, a significant heating of both the electronic and the lattice subsystem has been observed [18]. Motivated also by these evidences, we have recently addressed the issue of developing a theoretical model suitable for the microscopic simulation of realistic devices [19, 20], proposing an MC-based global kinetic approach [21]. In particular, our analysis has allowed to include on equal footing the fully three-dimensional nature of the multiband electronic structure as well as of the phonon degrees of freedom and is therefore able to overcome significant limitations affecting some of the MC techniques previously cited.

In the present chapter, we are presenting a review of the above-mentioned global MC approach. In particular, in Section 2 we shall describe the general coupled carrier-phonon model and the kinetic approach based on the MC sampling of the closed set of dynamical equations for the corresponding distribution functions; in Section 3, we will apply our model and method to a prototypical resonant-phonon THz QCL design, highlighting the interplay between the nonequilibrium electron and phonon populations; finally, in Section 4, we shall summarize and draw a few concluding remarks.

2. Physical model and kinetic description

First of all, to provide quantitative and predictive insight into the details of the carrier dynamics in unipolar multisubband devices such as QCLs, a proper treatment of the fully three-dimensional nature of the problem is mandatory [22]. In most configurations of practical interest, the description of the coupled electron-phonon nonequilibrium dynamics may safely be limited to a semiclassical approach for both subsystems. Within this scheme, the kinetic variables of interest are the electron and phonon distribution functions. The latter obey a coupled set of nonlinear equations of motion which may be formally derived in the framework of the quantum-kinetic theory, following a conventional single-particle correlation expansion procedure and neglecting all nondiagonal carrier and quasiparticle density matrix elements [23]. Coherent phonons and electron phase-coherence effects, which may show up within the transient ultrafast timescale [22], are indeed expected to play a minor role in the device’s steady-state operation regime considered here. Even though the system under investigation is spatially inhomogeneous, the LO phonon subsystem may be described in terms of a density matrix diagonal in q, parametrized by a meso/macroscopic spatial coordinate related to the thermal-transport space-scale. In other words, it could still be a good approximation to assume the phonon modes and interaction potentials of the host bulk material.

In principle, the proper treatment of the steady-state properties of our prototypical device should therefore be based on the set of coupled electron and phonon Boltzmann equations. This task is quite demanding, since it requires, for the latter, the inclusion of both acoustic and optical modes, ideally with finite size and quantization effects [12]. A simpler, though reasonable, starting point suggests to focus on the main physical aspects of the energy redistribution between charge carriers and lattice degrees of freedom. The basic idea is then to consider the full electron subsystem (i.e., the complete set of active region and injector subbands within each period) coupled via Frölich interaction with bulk LO phonon modes, and to include the decay of the latter into acoustic modes (lattice thermal bath at temperature $T_L$ ) via a phenomenological lifetime $\tau$ in their dynamical equation. The role of electron-acoustic phonon scattering is expected to be minor, due to the much lower rates, and for this reason it is not included in the present model.

In other words, energy dissipation is modeled—and basically occurs—as a two-step process: first, the carrier subsystem transfers a significant amount of energy to the LO phonon one; then, the latter transfers its excess energy to additional degrees of freedom (acoustic phonons) that are not taken into account dynamically and are characterized by a much larger heat capacity. This second step is described via a standard relaxation-time approximation. The structure of the closed set of coupled equations corresponding to the above-described dynamics is then the following:

$$\frac{\mathrm{d}{f}_{\nu {\mathbf{k}}_{\parallel }}}{\mathrm{d}t}={\left.\frac{\mathrm{d}{f}_{\nu {\mathbf{k}}_{\parallel }}}{\mathrm{d}t}\right|}_{\mathrm{e}-\mathrm{LO}}+{\left.\frac{\mathrm{d}{f}_{\nu {\mathbf{k}}_{\parallel }}}{\mathrm{d}t}\right|}_{\mathrm{s}}$$ E1

$$\frac{\mathrm{d}{n}_{\mathbf{q}}}{\mathrm{d}t}={\left.\frac{\mathrm{d}{n}_{\mathbf{q}}}{\mathrm{d}t}\right|}_{\mathrm{e}-\mathrm{LO}}+{\left.\frac{\mathrm{d}{n}_{\mathbf{q}}}{\mathrm{d}t}\right|}_{\text{loss}},$$ E2

where $f_{\nu {\mathbf{k}}_{\parallel }}$ is the carrier distribution function corresponding to the single-particle state in subband $\nu$ and with in-plane wavevector ${\mathbf{k}}_{\parallel }$ , while n_q is the average phonon occupation number of a single LO phonon mode with wavevector q.

Starting with the first term on the right-hand side of Eq. (1), in the framework of the Fermi’s golden rule approximation, the electron-phonon (e-LO) coupling hamiltonian produces the typical Boltzmann scattering structure. In particular, the latter consists of the following in- and out-scattering terms, coming from LO phonon emission and absorption processes:

$${\frac{\mathrm{d}}{\mathrm{d}t} f_{\nu {\mathbf{k}}_{\parallel }}}_{{|}_{\mathrm{e}}}=\sum _{\mathbf{q}\pm }\left(n_{\mathbf{q}}+\frac{1}{2}\pm \frac{1}{2}\right)\sum _{{\nu }^{\prime }{\mathbf{k}}_{\Vert }^{\prime }}\left[\left(1-f_{\nu {\mathbf{k}}_{\parallel }}\right)P_{\nu {\mathbf{k}}_{\parallel },{\nu }^{\prime }\mathbf{k}\text{'}}^{\mathbf{q}\pm }{f}_{{\nu }^{\prime }{\mathbf{k}}_{\Vert }^{\prime }}-\left(1-f_{{\nu }^{\prime }{\mathbf{k}}_{\Vert }^{\prime }}\right)P_{{\nu }^{\prime }{\mathbf{k}}_{\Vert }^{\prime },\nu \mathbf{k}}^{q\pm }{f}_{\nu {\mathbf{k}}_{\parallel }}\right].$$ E3

Here

$$P_{\nu {\mathbf{k}}_{\parallel },{\nu }^{\prime }{\mathbf{k}}_{\Vert }^{\prime }}^{\mathbf{q}\pm }=\frac{2\pi }{\mathrm{\hslash }}{\left|g_{\nu {\mathbf{k}}_{\parallel },{\nu }^{\prime }{\mathbf{k}}_{\parallel }^{\prime };\mathbf{q}}^{\pm }\right|}^2\delta \left({\epsilon }_{\nu {\mathbf{k}}_{\parallel }}-{\epsilon }_{{\nu }^{\prime }{\mathbf{k}}_{\parallel }^{\prime }}\pm {\epsilon }_{\mathbf{q}}\right),$$ E4

where $g_{\nu {\mathbf{k}}_{\parallel },{\nu }^{\prime }{\mathbf{k}}_{\parallel }^{\prime };\mathbf{q}}^{\pm }$ are the matrix elements of the q-dependent carrier-LO phonon coupling coefficients within the considered multisubband electronic basis $\nu {\mathbf{k}}_{\parallel }$ ; the sign + (−) refers to phonon emission (absorption) processes.

The second term, labeled with s, at the right-hand side of Eq. (1), generally accounts for further scattering mechanisms, of both intrinsic and extrinsic types, that may affect the carrier dynamics. Within the first sort of processes, for the kind of devices we are interested in, the electron-electron interaction has proven to generally play a significant role [3] and should be included in any realistic analysis (as will be done for the results presented in the next section). Regarding extrinsic mechanisms, such as, for example, interface roughness and electron-impurity scattering, they should be considered if one is interested in analyzing the behavior of a specific device. Indeed, such processes do not significantly modify the trend of the current-voltage characteristics since they have, unlike carrier-LO phonon scattering, a threshold-less nature and do poorly depend on the applied bias. On the contrary, they are strongly device/sample-dependent and their implementation therefore inevitably requires a phenomenological treatment. For these reasons, they will not be considered in the following.

The structure of Eq. (1) would allow for a self-consistent charge-conserving MC electron transport simulation on the condition that the phonon distribution n_q was known [22]. This typically occurs when the latter may be reasonably well approximated by the q-independent Bose-Einstein term corresponding to a given (quasi)equilibrium temperature. However, this is not the case we are presently interested in since the interplay between electron and phonon dynamics is here expected to significantly drive the LO phonon distribution out of equilibrium while the device is in operation. When nonequilibrium phonons are considered, their occupation numbers n_q are no longer q independent and have to be obtained by solving the corresponding dynamical equation, Eq. (2).

The Boltzmann transport equation for the electron subsystem, Eq. (1), with the previously described contribution in Eq. (3), goes then together with the phonon counterpart in Eq. (2). The two quantities at the right-hand side of the latter can be explicitly written as

$${\left.\frac{\mathrm{d}{n}_{\mathbf{q}}}{\mathrm{d}t}\right|}_{\mathrm{e}‐\mathrm{LO}}=\sum _{\pm }\pm \left(n_{\mathbf{q}}+\frac{1}{2}\pm \frac{1}{2}\right)\sum _{\nu {\mathbf{k}}_{\parallel },{\nu }^{\prime }{\mathbf{k}}_{\parallel }^{\prime }}\left(1-f_{\nu {\mathbf{k}}_{\parallel }}\right)P_{\nu {\mathbf{k}}_{\parallel },{\nu }^{\prime }\mathbf{k}\text{'}}^{\mathbf{q}\pm }{f}_{{\nu }^{\prime }{\mathbf{k}}_{\parallel }^{\prime }},$$ E5

$${\left.\frac{\mathrm{d}{n}_{\mathbf{q}}}{\mathrm{d}t}\right|}_{\text{loss}}=-\frac{n_{\mathbf{q}}-n_{\mathbf{q}}^{\mathrm{th}}}{\tau },$$ E6

with the $-\frac{n_{\mathbf{q}}-n_{\mathbf{q}}^{\mathrm{th}}}{\tau }$ term accounting for loss processes from the LO phonon subsystem due to anharmonic decay into acoustic modes, $n_{\mathbf{q}}^{\mathrm{th}}$ being the Bose-Einstein distribution at the lattice (acoustic) bath temperature $T_L$ .

The structure of the coupled set of sections, Eqs. (1) and (2), allows us to simulate both the electron and phonon dynamics by means of a MC particle-like technique. In particular, the latter consists of a generalized simulation approach suitable to describe a fixed number of electrons and a variable number of phonons, and can therefore be implemented in state-of-the-art MC simulation tools [22, 24, 25]. Technically, the q dependency of the phonon occupation numbers n_q can be managed by means of a combined self-scattering and rejection technique. In addition, the previously discussed symmetry-by-design of the device core region allows us to evaluate the device performances within a “closed-circuit” picture. More specifically, a periodic boundary condition scheme may be adopted in which the simulation/solution for the electron dynamics is limited to a single stage only. As described in detail in Ref. [24], every time an electron undergoes an inter-stage scattering process, it is properly re-injected into the simulated region and the corresponding charge contributes to the current through the device. This charge-conserving scheme allows for a purely kinetic evaluation of the device performances such as the gain spectrum or the current-voltage characteristics. The current density across the whole structure, for example, results as a pure output of the simulation just by a proper counting of the in- and out-stage scattering processes.

Within the above-described simulation framework, the LO phonon lifetime, $\tau$ , and the lattice temperature, $T_L$ , are the only free parameters. Realistic values for the former in bulk materials are in the range 6–9 ps [26], while the latter may be accessed by means of state-of-the-art microprobe band-to-band photoluminescence experiments, [17, 18].

3. Application to a prototypical THz coherent light source

As anticipated in the previous section, heating is a serious aspect generally affecting electronic and optoelectronic device performances at diverse stages. For the case of typical mid-infrared as well as THz resonant-phonon QCLs, a significant amount of LO phonons is generated in the core region during operation and a considerable amount of energy is dissipated in the device via Joule effect. Moreover, the combination of high LO phonon generation rates and limited thermal conductivity may cause a significant—and, in principle, detrimental—feedback on the electron relaxation kinetics and therefore on the device performances.

In recent years, experimental evidences highlighted the need for a deeper and more quantitative insight into the impact of hot-phonon phenomena on the nonequilibrium electron dynamics in THz QCL emitters [19, 20]. Motivated by these reasons, we have developed the global kinetic approach described in Section 2 [21]. In particular, the resonant-phonon device proposed in Ref. [27] appeared to be an interesting candidate deserving further investigation that will be reviewed here. As a starting point, besides the evident carrier-LO phonon scattering, also carrier-carrier interaction should be included in the electron dynamical equation, Eq. (1). This can be done employing the well-established time-dependent static-screening model which is commonly adopted in two-dimensional systems [28]. Regarding Eq. (2), the lifetime $\tau$ appearing there is set to the constant value of 6 ps throughout the whole simulated device operation range. This is a realistic guess based on the lower bound value in bulk materials [26] that allows us to investigate the relevance of nonequilibrium LO phonon effects already in the less favorable situation corresponding to their fastest decay into thermalized acoustic modes. Concerning the latter, the values for the input temperature parameter, $T_L$ , may be directly extracted from microprobe band-to-band photoluminescence analysis on the operating device. In particular, the experimental data evidence a significant device heating, with $T_L$ being bias dependent and locally higher than the cryostat temperature (100–180 K—varying with the applied bias—vs. 80 K) [19].

The global and three-dimensional kinetic approach reviewed here has the key feature of directly accessing the distribution functions of the electrons in the various subbands: these data are a straightforward output of the MC simulation. In other words, there is no a priori assumption of carrier intrasubband thermalization. The latter, instead, may or may not eventually show up because of the complex and sinergistic interplay between two-particle (carrier-carrier) and single-particle (carrier-LO phonon) interactions. It turns out that for the device that we are considering here [21], the nonequilibrium electron distributions do show a typical heated Maxwellian form. The corresponding subband effective temperatures are higher than the lattice ones, in good agreement with experimental findings [19]. As will be discussed later, a quite different scenario appears on the side of the LO phonon subsystem.

A key question when investigating and exploring novel QCL designs concerns the feedback that the nonequilibrium phonon dynamics has on the device gain performances. Since the electron cascading dynamics is strategically tailored and controlled by LO phonon emission processes, a gain reduction could be a reasonable expectation. However, any reliable/quantitative answer does require all the power and flexibility of our MC-based approach. In fact, a peculiar feature of MC simulations is that of allowing for truly simulated experiments, in which the diverse contributions to the dynamical equations may be selectively switched on and off, thus highlighting both their individual and synergistic roles. For the case of the prototypical THz QCL device considered here, when hot phonon effects are taken into account we indeed found a population inversion amounting to approximately 85% of the value corresponding to the ideal case in which the carrier subsystem interacts with a thermalized LO phonon population at the cryostat temperature [19]. A closer look at the subband population dynamics shows that this is in fact a “thermal backfilling” effect. More specifically, there is a significant reduction of the extraction efficiency out of the lower laser subband into the ground injector one, while the thermally activated depopulation of the upper laser states remains marginal. The presence of a nonequilibrium phonon population translates into an enhancement of the phonon absorption rates and this reduces the net energy relaxation/loss for the electrons [29].

Figure 1 shows the simulated current-voltage characteristics of the prototypical device we are considering. Indeed, results obtained from diverse MC simulated experiment—at diverse levels of description—are compared. The first set of data (triangles) corresponds to the plainest case, that is the one in which the complete (acoustic and optical) phonon subsystem is at equilibrium and behaves as a thermal bath at the cryostat (80 K) temperature. The other two datasets describe the electron transport performance when LO phonon nonequilibrium phenomena are dynamically taken into account and the acoustic phonon bath is at equilibrium either at the cryostat temperature (squares) or at higher (experimentally measured) ones (discs). All the three cases show the typical QCL behavior: as long as the injector and upper laser subbands remain aligned, the current density increases on increasing the bias; a negative differential resistance region then shows up at higher fields. When comparing the more realistic results (discs) with the other two, one finds out that the above-mentioned reduction in the electron cooling process out of the lower laser subband has a beneficial effect on the global transport dynamics across the device core region. Thermal backfilling actually prevents the accumulation of charge carriers in the injector ground subband, which is less efficiently coupled to the upper laser one, and forces electrons to follow alternative relaxation paths. Indeed, the resulting current increment may amount to more than 60% with respect to the equilibrium LO phonon case (triangles) and this is again a hot-phonon effect in the broader sense. Nonequilibrium LO phonon phenomena at the bare cryostat temperature for the acoustic phonon bath (squares) actually only account for a small fraction of the above-discussed findings.

Figure 1.

Simulated current density versus applied field characteristics, obtained under diverse conditions for the LO phonon subsystem and the acoustic phonon bath temperature. Triangles: LO and acoustic phonons at equilibrium at the cryostat (80 K) temperature; squares: nonequilibrium LO phonons and cryostat temperature for the acoustic phonon bath; discs: nonequilibrium LO phonons and measured $T_L$ for the acoustic phonon bath. Lines are a guide to the eye. Reprinted from [21].

A closer inspection of the current-voltage profiles allows for further remarks. In particular, the quantitative impact of the nonequilibrium phonon dynamics occurs at diverse levels, the most relevant being the specific device design and the particular operation condition, that is, the applied bias. From the theoretical modeling point of view, it follows that it is not possible to assume an a priori heated-phonon approximation parametrized by an effective temperature: the latter will unavoidably be device and bias dependent, and therefore fully equivalent to an a posteriori fitting parameter with no predictive added value.

The current density across the device, being a global quantity, somehow averages over the phonon distribution. In this respect, it does not unambiguously allow to discriminate between the effects of a quasiequilibrium population (thermal, though heated) and a nonthermal one. Indeed, the assessment of a nonequilibrium LO phonon population can and should be acknowledged by looking at the solution of Eq. (2). Figure 2 shows the n_q values obtained from the MC simulation, as a function of the in-plane modulus $q=\mid {\mathbf{q}}_{\text{in plane}}\mid$ . In particular, the plotted quantity is obtained averaging over the dependence on the out-of-plane wavevector component within the superlattice Brillouin minizone. These data (solid line) are compared with the thermal, q-independent, Bose distribution at the temperature $T_L$ (dashed line). Actually, the strictly quantitative values for the former and the latter depend on the applied bias; the qualitative behaviors at diverse bias show, however, a common feature: the build-up of a significant amount of small-q LO phonons within the core region of the operating device. This is indeed a characteristic fingerprint of the polar phonon emission process: the latter, whose rate decays as $1/q^2$ , privileges, at resonance, q wavevectors in the proximity of the zone center.

Figure 2.

Longitudinal optical-phonon occupation number (solid line) as a function of the in-plane wavevector modulus, for the operating THz QCL device considered in the text. The thermal Bose-Einstein distribution at the bias-dependent temperature $T_L$ is also shown for comparison (dashed line). Reprinted from [21].

The above-described nonthermal phonon population scenario has been recently probed and confirmed by microRaman spectroscopy studies. In particular, an extremely good agreement between simulated and measured data has been observed [20].

4. Summary and conclusions

In state-of-the-art nanomaterials and related optoelectronic nanodevices, the excited (i.e., nonequilibrium) charge carrier dynamics often involves the emission of a significant amount of phonons. As a consequence, the population of the latter may therefore be accordingly driven out of equilibrium as well. This is distinctly noticeable in new-generation device designs where relaxation via LO phonon emission is deliberately exploited as the key mechanism to achieve and maintain the desired transport/gain performances. In this contribution, we have analyzed hot electron versus hot phonon effects in a prototypical THz resonant-phonon QCL structure. In particular, we have reviewed a recently proposed MC-based global kinetic approach describing the entire interacting electron 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 modeled by a phenomenological lifetime. Energy dissipation occurs then via a two-step process: first, the strongly biased electron subsystem conveys a conspicuous amount of energy into the LO phonon one; then, the latter transfers the excess energy to additional degrees of freedom, having a much larger heat capacity. The results of our MC simulations are in very good agreement with measured data and highlight how and to what extent the nonequilibrium phonon dynamics may affect the electrooptical device performances.

The simulation scheme reviewed here may be extended to account for a more refined modeling of the LO phonon subsystem as well as to include a more realistic description of their decay into acoustic phonon modes and of the heat transport dynamics. Indeed, it is worth noting that the analysis presented so far, in principle, quantitatively applies only to a specific part of the device, the one whose $T_L$ temperature is selectively accessed and sampled by means of appropriate experimental techniques. Actually, mainly due to their multiinterface structure, the cross-plane thermal conductivity of these devices is greatly reduced with respect to the bulk materials, resulting in a nonuniform heating along the growth direction. The measured $T_L$ values are therefore local quantities that may show significant variations in the direction of the applied bias. A possible way to include in the proposed model these thermal resistance effects consists in assuming the value of $T_L$ —which is an input parameter—to be position dependent and set it from an experimental device temperature mapping.