Simulation of Electronic Quantum Devices: Failure of Semiclassical Models

: To simplify the design and optimization of new-generation nanomaterials and related electronic and optoelectronic quantum devices, energy dissipation versus decoherence phenomena are often simulated via local models based on the Wigner-function formalism. Such a local description is, however, intrinsically incompatible with the fully quantum-mechanical (i.e., non-local) nature of the dissipation-free carrier dynamics. While the limitations of such hybrid treatments have already been pointed out in the past in diverse contexts, the spirit of the present work is to provide a more cohesive and critical review. To this aim, we focus on the fundamental link between the Wigner-function picture and the density-matrix formalism. In particular, we show that, starting from well-established density-matrix-based models, the resulting Wigner-function dissipation and/or thermalization dynamics is necessarily non-local. This leads to the conclusion that the use of local Wigner function models borrowed from the semiclassical Boltzmann theory is formally not justiﬁed and may produce unreliable results, and that such simpliﬁed local treatments should be replaced by fully non-local quantum models derived, e.g., via the density-matrix formalism.


Introduction
Despite its intrinsically quantum-mechanical nature, the steady-state electro-optical response of a conduction electron within a bulk crystal may be safely described via semiclassical models based on the Boltzmann transport theory [1]. This fact has allowed the development of a wide variety of semiclassical simulation strategies [2], and, during the last decades, the latter have provided a strong impetus to the development of semiconductor physics and technology [3].
The replacement of semiclassical simulation schemes with genuine quantum-mechanical treatments implies, in general, a considerable increase in computational effort and resources; indeed, a microscopic treatment of various scattering mechanisms via proper Monte Carlo simulations [2] is often computationally too demanding already within the conventional Boltzmann theory; for quantum-transport simulations, the situation is even worse.
The design and optimization of new-generation nanomaterials as well as of related electronic and optoelectronic quantum devices is therefore frequently performed via partially phenomenological where the latter, in turn is the sum of a dissipation (diss) and of a thermalization (th) term, Here, the distance between the interfaces is l and z is the longitudinal transport direction. Reprinted from [67].
In the presence of slowly-varying electromagnetic forces and relatively large device sizes (compared to the typical electron de Broglie wavelength) the deterministic carrier dynamics can be safely described in classical terms via the conventional drift-diffusion model of the Boltzmann transport theory [2,18,20]. More specifically, adopting the customary parabolic-band approximation (with effective mass m * ) and in the presence of a static electric field only (i.e., no magnetic fields) described by a force F, the deterministic contribution in Equation (1) is given by: where v(k) =¯h k m * denotes the carrier group velocity. As anticipated, in a conventional semiconductor device, the deterministic carrier dynamics in Equation (3) is significantly modified by a variety of stochastic processes described by Equation (2); the latter, in turn, take place both within the device active region due to different scattering mechanisms (carrier-impurity, carrier-phonon, carrier-carrier, etc.) and close to the device spatial boundaries due to carrier exchange between the device and one or more external charge reservoirs. The explicit form of the corresponding dissipation and thermalization contributions in Equation (2) strongly depends on the nature of the transport model adopted. As we now discuss, in any semiclassical simulation scheme, dissipation as well as thermalization processes are always described via spatially local models, and the same applies to the Wigner-function modeling as well (see Section 3). Such local treatments can, however, be subdivided into two different classes: phenomenological models based on relaxation-time schemes where the total electronic charge within the device is not necessarily conserved (open-system paradigm) and kinetic models based on Boltzmann-type schemes where the total number of simulated carriers is always conserved (closed-system paradigm).

Phenomenological Dissipation and Thermalization Models
As far as energy-dissipation within the device active region is concerned, the simplest choice is the adoption of a phenomenological model based on the relaxation-time approximation (RTA) [18,20], i.e., where the relaxation of the carrier distribution in (r, k) toward its equilibrium value f • (r, k) is described in terms of a space-and momentum-dependent relaxation rate Γ(r, k) encoding all relevant scattering mechanisms characterizing the operational conditions of the device under investigation; the space and momentum/energy dependence of Γ(r, k) may be extracted from fully microscopic Monte Carlo simulations [2] or modeled via simplified Fermi's Golden rule treatments.
As far as the device-reservoir thermalization is concerned, within the semiclassical picture, the latter may be easily described via a net spatial separation between the device active region Ω and the external charge reservoirs (given by one or more electrical contacts/device terminals). Within this spirit, the simplest and physically more intuitive thermalization model is the so-called inflow-boundary-condition (IBC) scheme [2], where the carrier distribution function entering the device at its spatial boundary region coincides with the quasiequilibrium distribution in the external charge reservoirs. More specifically, denoting with r b the generic point of the device boundary surface σ b and with f b (r b , k) the carrier distribution function entering the device at the surface point r b , such a boundary condition can be easily fulfilled/implemented via the following thermalization operator [20] ∂ f (r, k) ∂t where is the incoming component of the carrier group velocity v(k) at the boundary point r b ( n(r b ) being the unit vector normal to the boundary surface σ b ). Indeed, according to its definition in Equation (6), v in is different from zero for incoming carriers only; moreover, in view of the Dirac delta function, the action of the thermalization operator in Equation (5) is limited to the device boundary surface σ b . As a result, the above thermalization scheme is able, as requested, to fix the incoming carrier distribution to that of the external reservoirs, i.e., f (r b , k in ) = f b (r b , k in ), and has no influence on the outgoing carrier flux. The simplest physically relevant implementation of the IBC scheme in Equation (5) is that of a one-dimensional system/device (r, k → z, k) of length l in thermal contact with a left and a right charge reservoir. As schematically shown in Figure 2, in this case, the boundary surface σ b of the three-dimensional case reduces to two points: z b = ±l/2. This amounts to fixing the value of the incoming distribution function f (z, k) at these two points, i.e., f (−l/2, k > 0) and f (+l/2, k < 0).
The conventional inflow or U boundary condition scheme adopted in semiclassical device modeling [20] for a one-dimensional problem. The value of the distribution function f (z, k) is specified at the boundaries z b (k) of the active region, i.e., f (−l/2, k > 0) and f (+l/2, k < 0) are fixed by the incoming/inflowing carrier distribution function. Reprinted from [67].
From a physical point of view, the action and role played by the IBC thermalization scheme in Equation (5) are significantly different from those of the dissipation term in Equation (4): while the latter modifies the carrier distribution in any point r inside the device volume Ω, the former simply fixes the value of the incoming distribution on the boundary surface σ b and has no effect inside the device. Despite these profound differences, the two models share a common formal structure, i.e., they are both RTA-type operators of the form: where, for the dissipation model in Equation (4), we have while, for the IBC thermalization model in Equation (5), we have It is worth stressing that, while for the dissipation model in Equation (8) the relaxation rate Γ acts as a simple algebraic multiplication, for the IBC thermalization model in Equation (9) the latter is an integral operator acting on our phase-space coordinates r, k.
As we discuss in Section 4, such formal similarity may help in understanding the common origin of diverse anomalous results obtained applying the local dissipation and thermalization models introduced thus far to quantum-transport simulations based on the Wigner-function formalism.

Kinetic Dissipation and Thermalization Models
To overcome the well-known limitations of the phenomenological dissipation model in Equation (4), the key step is to replace the RTA scheme with a Boltzmann collision term of the form [2,18,20] ∂ f (r, k) ∂t where P(r; k, k ) = (1 − f (r, k)) P 0 (r; k, k ) (11) denotes the low-density scattering rate P 0 in r for the generic transition k → k, weighted by the usual Pauli-blocking factor, and simply reduces to P 0 (r; k, k ) in the low-density limit ( f (r, k) → 0). The Boltzmann-type collision operator in Equation (10) provides a kinetic description of energy dissipation in terms of stochastic carrier transitions dictated by the low-density scattering rates P 0 (r; k, k ). Opposite to the phenomenological nature of the relaxation rate Γ(r, k) in Equation (4), such scattering rates can be microscopically derived in quantum-mechanical terms via the well-known Fermi's golden rule including all relevant interaction mechanisms. Moreover, such a stochastic dynamics may be efficiently and intuitively simulated via a wide variety of Monte Carlo approaches [2]. As anticipated, the Boltzmann-type kinetic treatment in Equation (10) is intrinsically charge conserving, which implies that in this case the possible time variation of the total carrier number inside the device active region Ω is due to device-reservoir thermalization processes only, customarily described via the phenomenological IBC scheme in Equation (5). As a result, most of the semiclassical device modeling is currently based on a kinetic treatment of dissipation phenomena inside the device properly combined with a phenomenological treatment of environment-induced carrier thermalization processes. This is typically accomplished via generalized (i.e., particle non-conserving) Monte Carlo simulations able to account for the presence of external reservoirs by means of corresponding injection versus loss terms [15].
As already stressed, the IBC thermalization model in Equation (5) is grounded on an artificial spatial separation between device and external contacts/reservoirs (see Figures 1 and 2); such an abrupt jump is, in general, not well justified, and one may expect that for some device operating conditions the distribution of incoming carriers may display nonequilibrium features as well. To properly describe the (continuous) device-reservoir transition, the key step is to include in the simulated region a significant portion of the electrical contacts as well; to this aim, in addition to the Boltzmann collision operator in Equation (10) (describing energy-dissipation phenomena inside the device), carrier thermalization due to the presence of the electrical contacts can be described again via a corresponding Boltzmann-type operator, whose thermalization rates are different from zero within the contact regions only. As a result, one deals with a kinetic treatment of both dissipation and thermalization phenomena, and the latter may be conveniently simulated via standard (i.e., particle conserving) Monte Carlo simulations [2]. Compared to the phenomenological IBC scheme, such kinetic treatments of the device-reservoir coupling are by far computationally more demanding. To overcome this limitation, an alternative Monte Carlo scheme, especially suited for the modeling of one-dimensional devices, has been also proposed [68,69]; the latter provides a solution of a Boltzmann-like equation defined over the single-particle electronic states of the device, and is thus more efficient compared to conventional phase-space (r, k) simulations.

Partially Quantum-Mechanical Device Modeling Based on the Wigner-Function Formalism
The Wigner function is introduced as the Weyl-Wigner transform of the single-particle density matrix. Even though the latter is defined over the very same phase-space r, k of the semiclassical distribution, and that very often the notation f (r, k) is used for the Wigner function as well, this new quantum distribution, while still being a real quantity, is no longer positive-definite. Its time evolution is described by the so-called Wigner transport equation. The latter can still be written analogously to Equation (1); due to this similarity, Monte Carlo simulation techniques for its solution (so-called Wigner Monte Carlo schemes [32,35]) have been developed. However, despite this formal resemblance, the explicit form of the deterministic as well as of the stochastic Wigner-function terms are expected to differ from their semiclassical counterparts. The deterministic contribution is the quantum-mechanical generalization of the diffusion-plus-drift term in Equation (3) and can be conveniently formulated in terms of Moyal brackets [70]. Indeed, denoting with V(r) the scalar potential corresponding to the electric force F, one has where is the nonlocal Weyl-Wigner superoperator corresponding to the scalar-potential profile [20]. Comparing Equation (12) with its semiclassical counterpart in Equation (3), we notice that the diffusion term remains unchanged while the quantum-mechanical version of the drift term is still local in r but non-local in k.
The stochastic contribution can still be cast in terms of dissipation as well as thermalization phenomena, as in Equation (2). It is then a common practice to treat these processes within the same local models previously described, namely the RTA model in Equation (7) and the Boltzmann collision operator in Equation (10). However, within the Wigner-function picture, this local approach suffers of intrinsic limitations. The latter is highlighted in the remainder of this section, where selected examples focusing on the analysis of energy dissipation versus decoherence phenomena [61] (Section 3.1) and on the device-reservoir carrier thermalization dynamics [67] (Section 3.2) are briefly reviewed.

Failure of Local Dissipation Models
The incompatibility between the local dissipation models of the semiclassical theory and the (non-local) deterministic Wigner-function dynamics in Equation (12) is already evident in modeling a basic prototypical quantum-well (QW) device such as the one sketched in Figure 3. The effects of device spatial boundaries as well as the role of device-reservoir thermalization phenomena may be neglected assuming semi-infinite barriers and limiting our investigation to QW bound states only. Due to confinement along the growth (z) direction, the system three-dimensional (3D) electronic states exhibit the usual subband structure. However, within the prescription of the RTA dissipation model in Equation (4), we neglect the in-plane phase-space coordinates and adopt an effective one-dimensional (1D) description: The set of single-particle quantum numbers then reduces to the discrete index (n) of the confined 1D state only and (r, k) → (z, k). Figure 3 shows the three bound energy levels and squared envelope functions for the model case of a GaAs/(Al,Ga)As 12 nm-thick QW, with a conduction band offset of 0.3 eV. In the following, we resume the results of simulated experiments performed on this nanosystem in the low-temperature regime, thus assuming that the QW lowest energy level (n = 1) corresponds to the thermal-equilibrium condition. In particular, we analyze the spatial carrier-density profile, along the growth direction, corresponding to the 1D Wigner function f (z, k), namely and investigate its time evolution, when the system is initially prepared in the first excited level (n = 2), by solving the 1D version of the Wigner transport equation (Equation (1)) with the quantum-mechanical deterministic term in Equation (12) and the conventional RTA scheme in Equation (4). In Figure 4, the dissipation dynamics induced on our prototypical QW system by three different versions of the RTA approach are compared. Within the model labeled M1, we assume a space-and energy-independent rate Γ 0 ; this is often called phenomenological or macroscopic RTA approach. Within the model M2, we employ a space-dependent but energy-independent rate Γ(z), assumed to vanish within a given region − a 2 < z < + a 2 , and equal to Γ 0 outside; this may partially mimic the case of material-dependent semiclassical scattering rates, such as in the case of complex multi-layer quantum heterostructures. Finally, within the model M3, we consider a space-independent but momentum/energy-dependent rate Γ(k), vanishing for |k| < k min = √ 2m * min /h and equal to Γ 0 otherwise; this may describe the typical energy-threshold scenario of electron-optical phonon scattering in a variety of nanomaterials and nanodevices, e.g., new-generation quantum-cascade lasers [15].
The dissipation dynamics resulting from M1 [ Figure 4a] agrees with the physically valid scenario in which the charge distribution in the QW ground state n = 1 increases, at the expenses of the one of the initial state n = 2 (thin solid curve); after 1.6 ps (dash-dotted curve) the electronic system is already pretty close to its equilibrium distribution (solid curve). In the case of M2 [ Figure 4b] and M3 [ Figure 4c], the scenario is substantially different: besides a significant slowdown of the excited-level decay process, negative carrier distributions arise; this crucial difference with respect to M1 is an unambiguous fingerprint of unphysical electronic states.  (14) is plotted as a function of position at different times (0 ps: thin solid curve; 0.8 ps: dashed curve; and 1.6 ps: dash-dotted curve) for models: M1 (a); M2 (b); and M3 (c). Here, the following parameters are employed: 1/Γ 0 = 0.8 ps for M1, 1/Γ 0 = 0.4 ps and a = 5 nm for M2, and 1/Γ 0 = 0.5 ps and min = 40 meV for M3. The equilibrium carrier density is also shown for comparison (solid curve). Reprinted from [61].
A closer inspection of the dissipation dynamics may be provided by looking at the carrier populations f n of the QW states. The latter are given by the diagonal terms of the single-particle density matrix, which can be obtained from the 1D Wigner function via the inverse of the 1D Weyl-Wigner transform. In Section 4, the link between the density matrix and the Wigner function is more extensively addressed. Here, we just anticipate that the populations f n can be obtained applying to the 1D Wigner function f (z, k) the inverse of the 1D Weyl-Wigner transform W nn (z, k), i.e., Figure 5 shows the time evolution of the populations of the three QW bound states ( f 1 , f 2 , and f 3 ) as well as of the continuum, the latter being defined as N cont = N tot − ∑ 3 n=1 f n , with N tot denoting the total number of carriers. As expected, model M1 [ Figure 5a] shows a simple exponential depopulation of the initial level n = 2 into the ground state n = 1, while level n = 3 and the continuum are not affected by the dissipation dynamics.
Such a plausible scenario is, however, lost in models M2 [ Figure 5b] and M3 [ Figure 5c], where we do observe: (i) a slowdown of the n = 2 state depopulation; (ii) negative population values for states n = 1 and n = 3; and (iii) particularly for M3, a marked charge transfer into the continuum in spite of our zero-temperature analysis. The unavoidable conclusion is that the conventional RTA dissipation model in Equation (4) does not preserve the positivity of the carrier populations f n , and induces a fictitious interlevel coupling; the latter, in turn, may also artificially generate interlevel phase coherence.
The unphysical results of M2 and M3 are not related to the charge non-conserving nature of the RTA model and/or to the effective 1D modeling of our nanodevice. In the following, we attest that analogous problems also arise when: (i) a fully 3D treatment of the prototypical QW nanostructure is adopted; and (ii) the RTA dissipation scheme in Equation (4) is replaced by the (charge-conserving) Boltzmann collision term in Equation (10). Concerning the latter, it has the well-established in-minus out-scattering structure and may also be written as with P in (r, k; r , k ) = δ(r − r )P(r; k, k ) (17) and Both superoperators, P in and P out , are local in r; P out is local in k as well.
The rates entering the Boltzmann collision term in Equation (10) are evaluated via the conventional Fermi's golden rule using as noninteracting states standard 3D plane waves, instead of the nanostructure single-particle wavefunctions, and assuming as main dissipation mechanism the carrier interaction with a bulk GaAs LO-phonon system [2]. The quantities P(r; k, k ) are therefore space-independent. Within such 3D description, the Weyl-Wigner phase-space of the QW device is given by r ≡ r , z and k ≡ k , k z . The effective 1D carrier-density profile along the growth direction (z) in Equation (14) is then now replaced by and the effective 1D level population in Equation (15) by the following subband population Here, f k n is the population of the generic carrier state in subband n with in-plane wavevector k , obtained via the 3D version (n → k n) of the Weyl-Wigner transform in Equation (15). Figure 6 shows the low-temperature and low-density energy-dissipation dynamics for the GaAs/(Al,Ga)As prototypical QW nanostructure of Figure 3, resulting from the fully 3D Boltzmann bulk model just described. The spatial carrier density profile [ Figure 3a] is now only marginally affected by negative-value regions, if compared to the unphysical results obtained via the RTA dissipation model (Figures 4 and 5). However, analogously to the RTA case, a marked slowdown of the excited-level decay is evident: after 1.6 ps (dash-dotted curve), the charge distribution is still far from equilibrium (solid curve). Such anomalous scenario is confirmed by an insight into the various subband populations (Figure 3b), showing: (i) a marked slowdown out of the initial subband n = 2; (ii) negative values in subband n = 3; and (iii) despite our zero-temperature analysis, considerable values in the continuum. A further close-up in the analysis may be achieved repeating the simulated experiment reported in Figure 6 with a different initial condition, e.g., the QW low-temperature thermal-equilibrium state n = 1. In this case, a correct scattering superoperator would leave the electronic system unmodified. However, the results, reported in Figure 7, show that the Boltzmann collision term in Equation (10) drives the electronic system out of equilibrium, populates the continuum, and induces negative populations of the excited QW subbands. These outcomes confirm the intrinsic limitations of local treatments, which are known to be intrinsically incompatible with any quantum-mechanical framework.
Indeed, both the RTA results in Figures 4 and 5 and the Boltzmann-type ones in Figures 6 and 7 are based on bulk-like scattering models, where the relaxation rate Γ(r, k) in Equation (4) and the semiclassical scattering probability P(r; k, k ) in Equation (10) do not account for the QW electronic subband structure. The two terms in the Wigner transport equation (Equation (1)) are then intrinsically conflicting: the deterministic one accounts for the QW subband structure via its confinement potential V(r), while the scattering one does not. This contradiction, recently pointed out also within the density-matrix formalism [71], is responsible for the dissipation slowdown as well as for the wrong thermalization dynamics. Indeed, the steady-state solution of the bulk-like Boltzmann collision term in Equation (10) is just the Fermi-Dirac distribution. The latter is space-independent and has nothing to do with the thermal-equilibrium Wigner function f • (r, k) of the quantum device, which is space-dependent and significantly different from zero within the QW region only. The same as in Figure 6 but considering as initial state the thermal-equilibrium state n = 1. Reprinted from [61].
As a final remark, we notice that the unphysical features of the low-density Boltzmann-type simulations just reviewed come out to be further amplified in the presence of high carrier concentrations. As discussed in detail in [64], in high-density conditions, the Wigner function may still be negative, but also greater than unity. This implies that, for high carrier concentrations, the role played by the Pauli factor (1 − f (r, k)) in Equation (11) may deviate significantly with respect to the conventional semiclassical scenario (where it is always bounded between 0 and 1), leading to highly counter-intuitive effects such as a sort of "negative dissipation".

Failure of Local Thermalization Models
In the presence of a significant device-reservoir thermal coupling, the primary goal of any quantum simulation scheme is the solution of the steady-state Wigner transport equation in particular, two crucial issues need to be assessed: (i) Is the solution of Equation (21) always unique? (ii) When this is the case, is it always physically acceptable?
As we now discuss, the answers to these questions depend crucially on the choice of the thermalization model as well as on the presence/absence of dissipation phenomena within the device active region. In particular, we now briefly review the detailed analysis originally presented in [67]. The latter has shown intrinsic limitations of Wigner-function simulations based on the conventional IBC thermalization scheme in Equation (5), providing the following answers: (i) in the coherent limit, the solution of the steady-state Wigner equation (compatible with given spatial boundary conditions) is definitely not unique; and (ii) when dissipation phenomena are taken into account, the solution, although unique, may be physically unacceptable.
It is useful to start with a simplified 1D model, having the notable advantage to allow for an exact analytical treatment; in particular: (i) the deterministic dynamics is described via a delta-like potential profile; (ii) the dissipation dynamics is neglected (so-called coherent limit); and (iii) the device-reservoir thermalization is described via the semiclassical IBC scheme in Equation (5), whose 1D version is schematically depicted in Figure 2. Within such a simplified scheme, it is possible to demonstrate [67] that the solution is not unique. This is a consequence of the double degeneracy of the 1D electronic scattering states corresponding to our potential profile, which in turn is due to the time-reversal symmetry of the dissipation-free dynamics.
In the presence of energy dissipation phenomena, in contrast, the time-reversal symmetry is lost, and the solution of the steady-state Wigner equation is found to always be unique. However, such mathematical solution does not necessarily correspond to a physically acceptable Wigner function, namely to the Weyl-Wigner transform of a genuine single-particle density matrix (see Section 4).
To better focus on this intrinsic limitation, let us consider a realistic GaAs-based nanostructure characterized by a rectangular potential barrier and assume a monoenergetic carrier injection from the left spatial boundary only. Moreover, we describe energy dissipation via the RTA scheme in Equation (4), with a space-and energy-independent relaxation rate Γ 0 = 1/τ.
The resulting spatial carrier-density profiles corresponding to the coherent limit (solid curve) as well as to two different values of the relaxation time [τ = 5 ps (dashed curve) and τ = 50 fs (dash-dotted curve)] are reported in Figure 8. In the coherent limit, a strongly anomalous feature shows up: the density profile comes out to be rigorously symmetric, despite the assumed left-injection into a significantly high potential barrier. The presence of energy dissipation mechanisms (dashed and dash-dotted curve) induces some (apparently more physical) spatial asymmetry. However, for the parameters considered here, the barrier transmission coefficient is so small (of the order of 10 −3 ) that no charge should be present in the right portion of the device (see also Figure 10 in Section 4.1). Moreover, in addition to this, all three density profiles display negative-value regions. This outcome is qualitatively similar to the one reported in [37] for the case of a cosine-like potential, thus confirming the physical limitations of conventional Wigner-function simulations based on the semiclassical IBC scheme. Again, similar to the case of the local dissipation models previously discussed, the physical origin of these anomalies is that the deterministic and the thermalization term entering the steady-state Wigner equation (Equation (21)) are intrinsically incompatible. While the former accounts for the QW subband structure via the nanostructure potential V(r) [see Equation (12)], the same does not apply to the semiclassical IBC scheme [see Equation (5)], where the carriers injected into the device active region are implicitly treated as free particles (i.e., plane waves), regardless of the presence of the confinement potential.

Fully Quantum-Mechanical Device Modeling Based on the Density-Matrix Formalism
In view of the intrinsic limitations of local dissipation as well as thermalization models discussed thus far, it is imperative to replace such semiclassical treatments with fully quantum-mechanical schemes. Despite of the level of accuracy/approximation of the latter, intimately related to the nanosystem as well as to the specific phenomenon under investigation, two basic requirements are, however, always mandatory: (i) the positive-definite character of the quantum-mechanical state has to be preserved at any time; and (ii) in the absence of external electro-optical excitations as well as of nonequilibrium carrier injection, the solution of the steady-state Wigner equation (Equation (21)) should coincide with the thermal-equilibrium state.
It is in general a highly non-trivial task to check the conformity to the above requirements; to this aim, it is crucial to exploit the close relation between the Wigner function and the corresponding single-particle density matrix [11,15]. More specifically, adopting the very same notation employed in [67], and by compactly labeling with α the set of relevant quantum numbers for the single-particle electronic states of our semiconductor nanostructure, the density matrix ρ α 1 α 2 (whose diagonal terms of the density matrix describe the population of the generic single-particle state α, while the off-diagonal terms describe the quantum-mechanical phase coherence (or polarization) between states α 1 and α 2 ) can be written in terms of the Wigner function f (r, k) as where denotes the well-known Weyl-Wigner transform [70], and φ α (r) the real-space wavefunction of the electronic state α.
Applying the inverse Weyl-Wigner transform equation (Equation (22)) to the original Wigner transport equations (Equations (1) and (2)), in the absence of external electromagnetic excitations, the latter can be easily translated into the density-matrix equation where α is the energy of state α, and As far as the deterministic dynamics is concerned, the density-matrix term in Equation (25) is by far simpler compared to its Wigner-function version in Equation (12): the diagonal terms (α 1 = α 2 ) do not change in time, while the non-diagonal terms (α 1 = α 2 ) simply rotate with frequency ( 1 − 2 )/h. As far as the dissipation and thermalization terms in Equation (26) are concerned, the latter may be conveniently described via well-established density-matrix models able to fulfill the two basic requirements mentioned above [72].
Based on the density-matrix picture just recalled, in the remainder of this section, we discuss in a more formal way the intrinsic limitations of the semiclassical treatments, and, starting from physically reliable and robust density-matrix models, we derive their Wigner-function versions given by fully non-local scattering superoperators.

Non-Local Relaxation-Time Models
As discussed in Section 2, both dissipation and thermalization processes are routinely described via the generic RTA scheme in Equation (7). To better clarify the intrinsic limitations of such local model within a Wigner-function simulation, it is useful to analyze its density-matrix version. To this aim, applying to Equation (7) the inverse Weyl-Wigner transform equation (Equation (22)), one gets [61] where α 1 α 2 denotes the inverse Weyl-Wigner transform of the equilibrium or quasiequilibrium distribution function f (r, k) in Equation (7). It is now worth comparing Equation (27) (obtained by applying a semiclassical RTA model to the Wigner function) to the equation obtained applying the RTA scheme straightforwardly to the density-matrix evolution, i.e., where ρ α 1 α 2 = f α 1 δ α 1 α 2 denotes the (diagonal) equilibrium or quasiequilibrium density matrix (with respect to our single-particle basis α). Here, Γ α can be regarded as a sort of state-dependent effective inverse life-time, involving all relevant dissipation/thermalization mechanisms acting on the carrier in state α; indeed, in the semiclassical limit (ρ α 1 α 2 = f α 1 δ α 1 α 2 ), the density-matrix RTA model in Equation (29) reduces to the well-known injection-loss scheme [15] of conventional device modeling.
In addition to the population dynamics in (30), the density-matrix RTA model in (29) describes the dissipation/thermalization-induced decay of the interlevel polarization, known as decoherence process. As one can see, while the superoperator in Equation (29) is diagonal in the single-particle basis α, the superoperator Γ in Equation (27) is, in general, non-diagonal in such basis. Moreover, while it is possible to show [72] that the density-matrix RTA model in Equation (29) is always positivity preserving, the same does not apply, in general, to the non-diagonal RTA model in Equation (27) [and therefore to its Wigner-function version in Equation (7)], as confirmed by the negative carrier populations in Figures 4 and 5). The only exception is the case of a space-and energy-independent relaxation rate, Γ(r, k) = Γ 0 , as well as a state-independent rate, Γ α = Γ 0 , for which the two RTA models in Equations (27) and (29) coincide. This is indeed confirmed by the convincing results in Figure 4a. The density-matrix analysis presented so far suggests to replace the highly problematic semiclassical-like RTA scheme in Equation (7) with a fully quantum-mechanical Wigner-function model. To this aim, our strategy, originally proposed in [61], is to start from the physically sound density-matrix RTA scheme in Equation (29) and apply to it the following Weyl-Wigner transform: The resulting Wigner-function model is given by where Γ(r, k; r , k ) = ∑ is a fully nonlocal RTA superoperator expressed in terms of the relaxation rates Γ α of the density-matrix theory, and f (r, k) is the Weyl-Wigner transform of the equilibrium or quasiequilibrium density matrix ρ α 1 α 2 .
In striking contrast to the standard RTA model in Equation (7), the generalized version in Equation (32) intrinsically ensures the positivity of the spatial charge density as well as of the single-particle state populations; indeed: (i) the density-matrix RTA model in Equation (29) is known to preserve the positive-definite character of the density matrix ρ α 1 α 2 ; and (ii) the action of the Weyl-Wigner transform in Equation (31) does not alter such property.
At this point, it is crucial to recall (see Section 2) that the original local RTA model in Equation (7) may be employed to describe both dissipation and thermalization phenomena, and the same applies to the density-matrix based RTA model in Equation (29) as well. Indeed, as discussed in more detail in [72], the generic relaxation rate Γ α may describe energy dissipation inside the device toward an equilibrium-state population ( , or device-reservoir thermalization toward the quasiequilibrium state population induced by one or more external charge reservoirs (Γ α = Γ res α , f α = f res α ). It follows that also the new RTA model in Equation (32) may be employed for the description of dissipation as well as thermalization phenomena, and that in both cases the action of the new RTA superoperator is always non-local.
To illustrate the quality and physical soundness of the proposed non-local RTA scheme in Equation (32), we focus first on the study of energy dissipation phenomena. To this end, the new non-local RTA model has been applied once again to the prototypical QW nanosystem in Figure 3 adopting again the very same effective 1D model introduced in Section 3.1. In this case, the in-plane carrier dynamics is neglected and the confined states of our QW system as well as the corresponding dissipation rates are labeled by the integer number n only, i.e., Γ diss α ≡ Γ diss n . In particular, Figure 9  Let us finally apply the proposed non-local RTA scheme in Equation (32) to the investigation of device-reservoir thermalization. To this end, let us consider again the 1D problem of Figure 8, namely the steady-state quantum transport across a rectangular barrier due to a monoenergetic carrier injection from left. In this case, the single-particle states α corresponding to the 1D rectangular potential profile are conventional scattering states labeled by their 1D carrier wavevector k, and the device-reservoir relaxation rates depend on k as well, i.e., Γ res α ≡ Γ res k . To mimic the monoenergetic-injection scenario of Figure 8, let us assume a monoenergetic carrier distribution in the left reservoir, i.e., f res k ∝ δ(k − k 0 ), k 0 denoting the wavevector corresponding to the injection energy 0 . As a result, in this case, the device-reservoir coupling is described via a single relaxation rate Γ res 0 ≡ Γ res k 0 . The latter can be regarded as the inverse of the so-called device transit time of the semiclassical theory [15], i.e., the ratio between the device length l and the injection carrier group velocity v 0 = (hk 0 )/m * , which for our simulation parameters is of the order of 100 fs. The results obtained via the new non-local RTA model in Equation (32) are reported in Figure 10, where we have assumed again a monoenergetic thermalization with 0 = 50 meV, and we have employed a device-reservoir relaxation time 1/Γ res 0 = 100 fs. Opposite to the highly non-physical results of the local IBC scheme in Figure 8, here, the charge distributions are always non-negative, thus confirming the positive-definite character of the proposed non-local RTA model in Equation (32). Moreover, regardless of the dissipation level, the amount of charge in the right portion of the device is now definitely negligible, in total agreement with the fact that for our simulation parameters the transmission coefficient across the rectangular barrier is extremely small (of the order of 10 −3 ). In addition, while in the coherent limit (solid curve) as well as in the weak dissipation regime (dashed curve) one deals with the typical interference scenario between incoming and reflected wavefunctions, in the strong-dissipation regime (dash-dotted curve), such interference process is partially suppressed due to dissipation induced decoherence [20].  Figure 10. The same as in Figure 8 but replacing the semiclassical IBC scheme based on the local RTA model in Equation (7) with the fully non-local model in Equation (32). To this end, the monoenergetic left injection of the IBC scheme employed in Figure 8 is replaced by a thermal coupling of the device with a corresponding monoenergetic distribution in the left reservoir via a relaxation time 1/Γ res 0 = 100 fs. The coherent-limit result (solid curve) and the weak-dissipation one (dashed curve) are basically coincident; the dash-dotted curve refers to the strong-dissipation regime.
Comparing once again the present convincing results with the highly problematic ones in Section 3.2, we also notice that, while there, in Figure 8, the curves corresponding to the coherent-limit result (solid line) and the weak-dissipation one (dashed line) differ significantly, here, in Figure 10, they basically coincide (the two profiles are hardly distinguishable). Recalling that the so-called coherent limit (τ → ∞) is concretely fulfilled when the dissipation time τ is much larger than any other characteristic time, and that for our system the only other time-scale is the thermalization time 1/Γ res 0 (of the order of 100 fs), it is clear that in the weak-dissipation regime (τ = 5 ps) the dissipation time is about fifty times larger than the thermalization one; this tells us that the weak-dissipation result (dashed curve) should closely resemble the coherent-limit one (solid curve). Such a physically sound scenario is fully confirmed by the results shown in Figure 10. This is instead not the case for Figure 8, where, as discussed in more detail in [67], the anomalous difference between the weak-dissipation result (dashed curve) and the coherent-limit solution (solid curve) is due to the fact that the latter is not unique. The coherent-limit result shown in Figure 8 (solid curve) should be therefore totally disregarded, and the dissipation-free scenario is well reproduced by the weak-dissipation result (dashed curve).
To confirm the quality and physical soundness of the non-local RTA model in Equation (32), we have repeated the simulated experiments of Figure 10 for a thinner (w = 6 nm) and lower (V 0 = 75 meV) barrier. Opposite to the original (almost transmission-free) barrier profile in Figure 10, in this case, the transmission coefficient is about 25 %, leading to a significant tunneling dynamics across the barrier. This scenario is fully confirmed by the new results reported in Figure 11: while in the left portion of the device the role played by dissipation phenomena is similar to the one reported in Figure 10, here, in addition, we deal with a relevant tunneling across the barrier, leading to a significant charge distribution in the right portion of the device as well. In particular, while the coherent-limit result (solid curve) coincides again with the weak-dissipation one (dashed curve), in the strong-dissipation regime (dash-dotted curve), we deal again with a suppression of the left-side interference pattern, and to a decoherence-induced suppression of the tunneling dynamics.  Figure 11. The same as in Figure 10 but for a thinner (w = 6 nm) and lower (V 0 = 75 meV) barrier (see text).
It is finally important to stress once again that the physically sound results in Figures 9-11 obtained employing the new RTA model in Equation (32) are primarily due to its non-local character. To better elucidate the link between the non-local RTA model in Equation (32) and its local version in Equation (7), it is useful to realize that the former reduces to the latter in the presence of a fully local relaxation function, namely: Γ(r, k; r , k ) = Γ(r, k)δ(r − r )δ(k − k ) .
Recalling that, for the thermalization IBC scheme in Equation (5), we have Γ(r, k) = v in (r b , k)δ(r − r b ), the explicit form of the local relaxation function in Equation (34) is given by: This clearly shows that the local action of the conventional IBC scheme takes place on the device boundary region r b only, in clear contrast to the action of the non-local relaxation function in Equation (33). Indeed, regarding again any RTA model as an injection/source minus a loss/drain term [see Equation (30)], while in the conventional IBC scheme carrier injection versus loss takes place on the device spatial boundaries only, in the non-local model, the action of both injection and loss phenomena due to the presence of external charge reservoirs is present inside the device as well.

Non-Local Boltzmann-Type Models
The very same generalization procedure employed to derive the nonlocal RTA model in Equation (32) can also be adopted to overcome the limitations of the Boltzmann-like treatment pointed out in Section 3.1 (see Figures 6 and 7). Instead of defining the scattering superoperator directly within the Weyl-Wigner phase-space (r, k), one can start from a reliable dissipation model for the density matrix, and then translate it into the Wigner-function picture.
To this aim, we refer to the nonlinear density-matrix approach proposed in [73], which: (i) applies to any generic nanostructure; (ii) accounts for high-density effects; (iii) provides the correct thermal-equilibrium state; and (iv) preserves the positive-definite character of the single-particle density matrix. More specifically, for both carrier-phonon and carrier-carrier interaction mechanisms, energy dissipation and decoherence is described by the following nonlinear scattering superoperator where P α 1 α 2 ,α 1 α 2 are generalized scattering rates, whose explicit form is given in [73]. The nonlinearity factors (δ α 1 α 2 − ρ α 1 α 2 ) can be regarded as the quantum-mechanical generalization of the Pauli factors of the conventional Boltzmann theory [see Equation (10)].
The key step is then to apply to the density-matrix scattering superoperator in Equation (36) the Weyl-Wigner transform in Equation (31). The resulting Wigner-function scattering superoperator still has the in-minus-out structure of Equation (16), with the local terms of the semiclassical theory in Equations (17) and (18) replaced with the following nonlocal generalizations.
and P out (r , k ; r, k; r , It is worth stressing that, despite the strong formal similarity with the conventional Boltzmann transport theory, the generalized Wigner-function scattering rates in Equation (37) are not necessarily positive-definite.
The proposed quantum-mechanical generalization of the standard Boltzmann collision term in Equation (10) is thus intrinsically nonlocal. In particular, comparing Equation (37) with its semiclassical counterpart in Equation (11), it is evident that the action of the Pauli exclusion principle within the Wigner phase-space is itself nonlocal: the generalized in and out scattering rates for a given transition r, k → r , k depend on the value of the Wigner function in any other phase-space point r , k via the Pauli factor 1 − f (r , k ). A detailed investigation of such Pauli-blocking nonlocality is outside the scope of the present work and can be found in [64].
In the low-density limit f (r, k) → 0, the proposed scattering model in Equation (37) then reduces to: and P out (r, k; r , k ) = 1 (2π) 3 ∑ α 1 α 2 α 1 α 2 W α 1 α 2 (r, k)P * To test the suitability of the proposed nonlocal Boltzmann superoperator, we have applied the latter to the case of the simulated experiment of Figure 6; the resulting energy-dissipation dynamics is presented in Figure 12. Here, opposite to the unphysical behaviors obtained via the semiclassical Boltzmann model (see Figures 6 and 7), both the spatial carrier distributions (Figure 12a) and the corresponding subband populations (Figure 12b) are always positive-definite. Moreover, analogously to the physically sound results of the RTA model M1 [see Figures 4a and 5a], the system now correctly relaxes toward the equilibrium state (solid curve in panel a) without unphysical interlevel couplings and continuum populations. Here, due to the trace-preserving character of the original density-matrix model in Equation (36), the proposed nonlocal Boltzmann-type superoperator in Equation (37) is intrinsically charge-conserving. cont Figure 12. The same as in Figure 6 but replacing the local Boltzmann model in Equation (10) with the proposed nonlocal model in Equation (37). Reprinted from [61].
As pointed out in Section 2.2, a possible strategy to improve the quality of semiclassical device modeling is to replace the standard IBC thermalization scheme with a global simulation of the device-plus-reservoir system. To this aim, an efficient approach is to describe reservoir-induced carrier thermalization via an effective Boltzmann-type collision term, as originally proposed in [68,69]. It is worth stressing that the very same strategy can be adopted for quantum-simulation schemes as well; indeed, instead of accounting for thermalization phenomena via the non-local RTA scheme in Equation (32), it is possible, in principle, to describe reservoir-induced carrier thermalization via the non-local Boltzmann-type superoperator in Equation (16) equipped with corresponding device-reservoir thermalization rates formally identical to the non-local dissipation rates in Equation (37).

Summary and Conclusions
The simulation of novel nanomaterials and of related electronic quantum devices is often grounded on partially phenomenological transport models based on the Wigner-function formalism. Within these simplified modeling strategies, energy dissipation versus decoherence phenomena are still described via local models borrowed from the semiclassical Boltzmann theory, namely relaxation-time treatments or Boltzmann-type scattering approaches. Such a semiclassical treatment is intrinsically incompatible with the fully quantum-mechanical coherent (i.e., dissipation-free) carrier dynamics within the device active region. Indeed, during the last decades, the intrinsic limitations of such hybrid treatments emerged in different contexts, both for semiclassical thermalization models based on the conventional boundary-condition scheme and for dissipation versus decoherence treatments based on the relaxation-time approximation, as well as on Boltzmann-like scattering models.
In this paper, our primary goal is twofold: (i) We review selected simulated experiments, showing the above-mentioned limitations in different contexts, and providing a cohesive theoretical treatment able to clearly identify the origin of such diverse anomalies faced within the Wigner-function picture.
(ii) Starting from the density-matrix formalism, we show how to derive genuine quantum mechanical dissipation models, able to overcome the intrinsic limitations just reviewed.
The main conclusion of our analysis is as follows: within the Wigner-function picture, a correct treatment of both dissipation versus decoherence phenomena inside the device and of the device-reservoir thermalization process requires the use of intrinsically non-local scattering models. In addition, our investigation shows that the most natural theoretical framework for a fully quantum-mechanical simulation of electronic and optoelectronic nanodevices is the density-matrix picture, which comes out to also be computationally less expensive. Indeed, in the presence of non-local dissipation and/or thermalization models, a direct numerical solution of the corresponding transport equation via finite-difference schemes over the Wigner phase-space is by far computationally more demanding (compared to the local case), and a direct numerical solution of the density-matrix equation comes out to be, in general, less expensive. No matter which is the adopted computational strategy (Wigner-function versus density-matrix solution), at any time, one can simply switch between the two pictures via the corresponding Weyl-Wigner transform.