Energy Dissipation and Decoherence in Solid-State Quantum Devices: Markovian versus non-Markovian Treatments

Abstract

The design and optimization of new-generation solid-state quantum hardware absolutely requires reliable dissipation versus decoherence models. Depending on the device operational condition, the latter may range from Markov-type schemes (both phenomenological- and microscopic- like) to quantum-kinetic approaches. The primary goal of this paper is to review in a cohesive way virtues versus limitations of the most popular approaches, focussing on a few critical issues recently pointed out (see, e.g., Phys. Rev. B 90, 125140 (2014); Eur. Phys. J. B 90, 250 (2017)) and linking them within a common framework. By means of properly designed simulated experiments of a prototypical quantum-dot nanostructure (described via a two-level electronic system coupled to a phonon bath), we shall show that both conventional (i.e., non-Lindblad) Markov models and density-matrix-based non-Markov approaches (i.e., quantum-kinetic treatments) may lead to significant positivity violations. While for the former case the problem is easily avoidable by choosing genuine Lindblad-type dissipation models, for the latter, a general strategy is still missing.

1. Introduction

Quantum-mechanical state superposition and correlation (i.e., entanglement) are the key players for the concrete realization of quantum information processing devices [1,2]. In particular, the pivotal ingredient for many solid-state implementations is electronic phase coherence [3,4]. The latter is, however, strongly hindered by dissipation and/or decoherence phenomena [5]. For the design and optimization of new-generation electronic quantum devices [6,7], it is then imperative to employ reliable dissipation versus decoherence models. The most popular quantum-mechanical picture for the description of phase coherence versus dissipation/decoherence in open quantum systems is the well-known density-matrix formalism [8], recalled and applied to the case of semiconductor nanodevices in Section 2. Within this picture, the simplest example is a two-level system, described by the following (two-by-two) density matrix:

$$\left(\begin{array}{cc}{\rho }_{bb}& {\rho }_{ba}\\ {\rho }_{ab}& {\rho }_{aa}\end{array}\right)=\left(\begin{array}{cc}{f}_b& p\\ p^{*}& f_a\end{array}\right)$$

(1)

in terms of the ground- and excited-level populations, $f_a$ and $f_b$ respectively, as well as of the interlevel phase coherence (or polarization) p. In this case, the simplest dissipation model is the well-known $T_1/T_2$ scheme [7]. This fully phenomenological approach accounts for dissipation versus decoherence processes via bare relaxation-time approximations: the evolution of the density matrix in Equation (1) is described by the following set of coupled equations:

$$\frac{d{f}_{a/b}}{dt}=-\frac{f_{a/b}-f_{a/b}^{°}}{T_1},\frac{dp}{dt}=\frac{{\Delta }_{\mathrm{c}}}{ıħ}p-\frac{p}{T_2}$$

(2)

(${\Delta }_{\mathrm{c}}={ϵ}_b-{ϵ}_a$ denoting the interlevel energy splitting), for which the solution is simply given by

$$f_{a/b}\left(t\right)=f_{a/b}^{°}+\left(f_{a/b}\left(0\right)\right)e^{-\frac{t}{T_1}},p\left(t\right)=p\left(0\right)e^{\frac{{\Delta }_{\mathrm{c}}t}{ıħ}}{e}^{-\frac{t}{T_2}}.$$

(3)

At long times, the level populations $f_{a/b}$ tend to their thermal-equilibrium values $f_{a/b}^{°}$ with a relaxation time $T_1$, while the interlevel polarization p decays to zero with another relaxation time, $T_2$. In spite of its success in the interpretation of many ultrafast optical experiments [9], the $T_1/T_2$ model may lead to totally unphysical results when the two relaxation times are treated as independent parameters.

Indeed, as discussed below, a crucial prerequisite of any reliable dissipation model is to preserve the positivity of the density matrix, namely the positive-definite character of its eigenvalues. For the case of our two-level system, such basic requirement translates into the condition that the eigenvalues

$${\mathsf{\Lambda }}_{\pm }=\frac{1}{2}\left((f_a+f_b)\pm \sqrt{{(f_a-f_b)}^2+4{\left|p\right|}^2}\right)$$

(4)

of the density matrix in Equation (1) are nonnegative or, equivalently, that the determinant $f_a{f}_b-{\left|p\right|}^2$ is nonnegative. Positivity violation therefore does not necessarily correspond to a negative value of either $f_a$ or $f_b$, as it may also occur in the presence of nonnegative level populations. This is the case, for example, of the $T_1/T_2$ solution in Equation (3) for particular combinations of the two phenomenological parameters, as shown in Figure 1. This pathological behavior is well known and is ascribed to the fact that the two relaxation times $T_1$ and $T_2$ are not physically independent quantities; indeed, by adopting more refined descriptions [7], they may be expressed in terms of microscopic scattering rates corresponding to the various interaction processes. The first obvious conclusion drawn from the above scenario is the potential inadequacy of barely phenomenological models and therefore the need for fully microscopic treatments.

The microscopic derivation of reliable scattering superoperators is one of the most challenging problems in quantum physics. For purely atomic and/or photonic systems, dissipation and decoherence phenomena may successfully be treated via adiabatic-decoupling schemes [8] in terms of extremely simplified models based on a few key parameters; within such effective treatments, the main goal is to derive a suitable form of the system-environment Liouville superoperator, able to preserve the positivity of the system density matrix [10]. This is typically accomplished by identifying so-called Lindblad superoperators [11] expressed in terms of a few crucial system–environment coupling constants. In contrast, for solid-state materials and related devices, the complex many-body quantum evolution results in a highly nontrivial interplay between electronic phase coherence and dissipation/decoherence [12,13,14,15,16], thus requiring the adoption of microscopic treatments.

Generally speaking, the microscopic derivation of solid-state dissipation models within the density-matrix picture may involve one or more of the following three key steps: (i) mean-field approximation, (ii) adiabatic or Markov limit, and (iii) semiclassical or diagonal approximation.

When all these three steps/approximations are performed, the collision term familiar from the Boltzmann theory is obtained; the latter, when applicable, constitutes a robust/reliable particle-like description in purely stochastic terms, thus providing in any case physically acceptable results.

In contrast, the combination of the first two steps only, namely mean-field approximation and adiabatic limit, allows one to derive so-called Markovian scattering superoperators, for which action may lead again to positivity violations [17]. Indeed, as originally pointed out by Spohn and coworkers [18], the choice of the adiabatic decoupling strategy is definitely not unique. Only the case studied by Davies [10], namely a “small” subsystem interacting with a thermal environment, could be shown to preserve positivity. However, such result was restricted to finite-dimensional subsystems (e.g., few-level atoms) and to the particular projection scheme of the partial trace. Thus, as such, it cannot be straightforwardly extended to the study of solid-state systems.

To overcome this severe limitation, a few years ago, an alternative and more general Markov procedure has been proposed [19]; the latter allows for a microscopic derivation of Lindblad-type scattering superoperators [11], thus preserving the positive-definite nature of the electronic density matrix. More recently, such an alternative Markov scheme combined with the conventional mean-field approximation has allowed for the derivation of a positive-definite nonlinear equation for the single-particle density matrix [20,21], able to describe both carrier–phonon and carrier–carrier interaction; the latter has been recently applied to the investigation of scattering nonlocality in GaN-based materials [22] and carbon nanotubes [23] as well as to the study of carrier capture processes [24,25].

For strong (i.e., non-perturbative) system–environment couplings combined with extremely short excitation and/or detection timescales, the adoption of the adiabatic or Markov schemes just recalled becomes questionable, and memory effects may be investigated within the density-matrix formalism via so-called quantum-kinetic schemes [12,15]. Indeed, stimulated by the pioneering papers by Haug and coworkers [26] as well as by Kuhn and coworkers [27], over the last decade, several groups have routinely employed such non-Markovian techniques to study a wide spectrum of ultrafast coherent phenomena in semiconductor bulk and nanostructures [28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61]. In spite of the undoubted success, these quantum-kinetic treatments, based on the mean-field approximation only, may lead once again to positivity violations. Indeed, such potential limitation, originally pointed out in the early days of electron-phonon quantum kinetics by Zimmermann and coworkers [62], has been recently investigated in more detail [63,64].

The primary goal of this paper is to review in a cohesive way virtues versus limitations of the most used dissipation quantum models employed in the simulation of state-of-the-art electronic and optoelectronic nanodevices, focussing on a few critical issues recently pointed out in Reference [21] as well as in Reference [64] and linking them within a common framework. More specifically, by means of properly designed simulated experiments of a prototypical quantum-dot nanostructure (described via an electronic two-level system coupled to a phonon bath), we shall show the following:

(i)

conventional (i.e., non-Lindblad) Markov models may lead to significant positivity violations;

(ii)

such intrinsic limitations may be avoided adopting properly designed Lindblad-type Markov schemes;

(iii)

density-matrix-based non-Markov models, namely quantum-kinetic treatments, may lead to positivity violations as well.

The most important conclusion of our investigation is that the presence of positivity violations is ascribed not only to the adiabatic approximation but also to the mean-field approximation. Indeed, while for the case of Markov treatments (based on both adiabatic and mean-field approximations), the problem is easily avoidable by employing genuine Lindblad-type dissipation models, for the case of non-Markov treatments (based on the mean-field approximation only), a general strategy is still missing. As a result, the unusual conclusion is that, in this case, two approximations work better than one.

The paper is organized as follows: In Section 2, we shall focus on Markovian models, comparing the conventional adiabatic-decoupling scheme with the Lindblad-type one and generalizing the latter to the nonlinear (i.e., degenerate) regime. In Section 3, we shall discuss virtues versus limitations of density-matrix-based non-Markovian models, namely quantum-kinetic treatments, pointing out specific conditions/regimes which may lead to positivity violations. Finally, in Section 4, we shall summarize and draw a few conclusions.

2. Markovian Dissipation Models

Within the spirit of the usual perturbation theory, the global solid-state Hamiltonian (electrons plus various crystal excitations, e.g., phonons, plasmons, etc.) may be schematically written as

$$\widehat{\mathbf{H}}={\widehat{\mathbf{H}}}_{°}+\sum _s{\widehat{\mathbf{H}}}_s^{\prime },$$

(5)

where the first term is the unperturbed contribution that can be treated exactly and the second term describes a number of perturbations ${\widehat{\mathbf{H}}}_s^{\prime }$, corresponding to various interaction mechanisms (e.g., carrier–phonon, carrier–carrier, etc.), which are typically treated within some approximation scheme.

2.1. Conventional Adiabatic-Decoupling Scheme

Following the fully operatorial approach originally proposed in Reference [17]) and described in more detail in Reference [7], the second-order (or incoherent) contribution to the time evolution of the global (e.g., carriers plus phonons) density-matrix operator $\widehat{\rho }$ obtained via the conventional adiabatic limit can be written as

$${\left.\frac{d\widehat{\rho }}{dt}\right|}_{\mathrm{inco}}=\frac{1}{2}\sum _s\left({\widehat{\mathbf{a}}}^s\widehat{\rho }{\widehat{\mathbf{b}}}^{s†}-{\widehat{\mathbf{a}}}^{s†}{\widehat{\mathbf{b}}}^s\widehat{\rho }\right)+\mathrm{H}.\mathrm{c}.,$$

(6)

where

$${\widehat{\mathbf{a}}}^s=\frac{{\widehat{\mathbf{H}}}_s^{\prime }}{ħ},{\widehat{\mathbf{b}}}^s=\frac{1}{ħ}{\int }_{-\infty }^{+\infty }{e}^{-\frac{{\widehat{\mathbf{H}}}_{°}{t}^{\prime }}{ıħ}}{\widehat{\mathbf{H}}}_s^{\prime }{e}^{\frac{{\widehat{\mathbf{H}}}_{°}{t}^{\prime }}{ıħ}}d{t}^{\prime },$$

(7)

and H.c. denotes the Hermitian conjugate.

The above Markov evolution is definitely non-Lindblad and therefore does not necessarily preserve the positivity of the global density matrix $\widehat{\rho }$ [19]. The scattering superoperator in Equation (6) may suitably be expressed in terms of generalized scattering rates. More specifically, denoting with $|i\rangle$ the generic eigenstate of the noninteracting Hamiltonian ${\widehat{\mathbf{H}}}_{°}$ and with ${ϵ}_i$ the corresponding energy level, one gets

$${\left.\frac{d{\rho }_{i_1{i}_2}}{dt}\right|}_{\mathrm{inco}}=\frac{1}{2}\sum _{s,i_1^{\prime }{i}_2^{\prime }}\left({\mathcal{P}}_{i_1{i}_2,i_1^{\prime }{i}_2^{\prime }}^s{\rho }_{i_1^{\prime }{i}_2^{\prime }}-{\mathcal{P}}_{i_1^{\prime }{i}_1^{\prime },i_1{i}_2^{\prime }}^{s*}{\rho }_{i_2^{\prime }{i}_2}\right)+\mathrm{H}.\mathrm{c}.$$

(8)

with generalized scattering rates

$${\mathcal{P}}_{i_1{i}_2,i_1^{\prime }{i}_2^{\prime }}^s={\mathbf{a}}_{i_1{i}_1^{\prime }}^s{\mathbf{b}}_{i_2{i}_2^{\prime }}^{s*}.$$

(9)

It is possible to show [7] that their diagonal (i.e., semiclassical) elements ($i_1{i}_1^{\prime }=i_2{i}_2^{\prime }$) coincide with the standard Fermi’s-golden-rule prescription:

$${\mathcal{P}}_{ii,i^{\prime }{i}^{\prime }}^s=\frac{2\pi }{ħ}{\left|\langle i\left|{\widehat{\mathbf{H}}}_s^{\prime }\right|i^{\prime }\rangle \right|}^2\delta ({ϵ}_i-{ϵ}_{i^{\prime }}).$$

(10)

Since the study of electro-optical processes in solid-state systems mainly relies on physical quantities that depend on the electronic-subsystem coordinates only, it is customary to introduce a many-body density-matrix operator

$${\widehat{\rho }}_{\mathrm{c}}=\mathrm{tr}{\left\{\widehat{\rho }\right\}}_{\mathrm{p}},$$

(11)

where the nonrelevant phononic (p) degrees of freedom have been traced out of the global density-matrix operator $\widehat{\rho }$. It is worth stressing that such treatment of carrier–phonon interaction applies to other bosonic degrees of freedom as well (e.g., photons, plasmons, etc.).

More specifically, by denoting with ${\widehat{\rho }}_{\mathrm{p}}^{°}$ the equilibrium density-matrix operator of the phononic subsystem and by adopting a carrier–phonon mean-field approximation via the following state factorization

$$\widehat{\rho }={\widehat{\rho }}_{\mathrm{c}}\otimes {\widehat{\rho }}_{\mathrm{p}}^{°},$$

(12)

it is possible to show [7] that the reduced dynamics dictated by the global evolution in Equation (6) is still of the same form:

$${\left.\frac{d{\widehat{\rho }}_{\mathrm{c}}}{dt}\right|}_{\mathrm{inco}}=\frac{1}{2}\sum _s\left({\widehat{\mathbf{a}}}_{\mathrm{c}}^s{\widehat{\rho }}_{\mathrm{c}}{\widehat{\mathbf{b}}}_{\mathrm{c}}^{s†}-{\widehat{\mathbf{a}}}_{\mathrm{c}}^{s†}{\widehat{\mathbf{b}}}_{\mathrm{c}}^s{\widehat{\rho }}_{\mathrm{c}}\right)+\mathrm{H}.\mathrm{c}..$$

(13)

Here, the explicit form of the reduced or electronic operators ${\widehat{\mathbf{a}}}_{\mathrm{c}}^s$ and ${\widehat{\mathbf{b}}}_{\mathrm{c}}^s$ can be derived starting from the global scattering operators ${\widehat{\mathbf{a}}}^s$ and ${\widehat{\mathbf{b}}}^s$ in Equation (7). We stress that, in spite of their very same formal structure, Equations (6) and (13) describe the system dynamics at different levels. This is confirmed by the fact that, while the global operators in Equation (7) are always Hermitian, the electronic ones, ${\widehat{\mathbf{a}}}_{\mathrm{c}}^s$ and ${\widehat{\mathbf{b}}}_{\mathrm{c}}^s$, are generally non-Hermitian, a clear fingerprint of dissipation-versus-decoherence processes induced by the phononic subsystem on the carrier one.

Within the above description, although a statistical average over the phononic degrees of freedom has been performed, the electronic subsystem is still treated via a many-body picture. Nevertheless, in the investigation of solid-state quantum materials and related devices, many of the physical quantities of interest are described via single-particle electronic operators. This suggests to treat the electronic subsystem via an additional mean-field (or Hartree–Fock) approximation. When applicable [7], this last approximation step allows one to get nonlinear single-particle dissipation models (see Section 2.3). In the low-density limit/regime, however, carrier–carrier interaction as well as Pauli-blocking effects can safely be neglected, and the many-body Hartree–Fock scheme just mentioned can conveniently be replaced by a simple one-electron model, i.e., the properties of a system of noninteracting electrons are fully described by one electron only. This amounts to treating the carrier subsystem in terms of a one-electron density matrix ${\rho }_{{\alpha }_1{\alpha }_2}$, i.e.,

$${\widehat{\rho }}_{\mathrm{c}}=\sum _{{\alpha }_1{\alpha }_2}|{\alpha }_1\rangle {\rho }_{{\alpha }_1{\alpha }_2}\langle {\alpha }_2|,$$

(14)

where $|\alpha \rangle$ denotes the generic one-electron eigenstate.

In the particular yet physically relevant case of a low-density carrier gas (c) interacting with a phonon bath (p), the explicit form of the global Hamiltonian in Equation (5) is given by

$$\widehat{\mathbf{H}}={\widehat{\mathbf{H}}}_{°}+{\widehat{\mathbf{H}}}^{\prime }=\left({\widehat{\mathbf{H}}}_{\mathrm{c}}+{\widehat{\mathbf{H}}}_{\mathrm{p}}\right)+\left({\widehat{\mathbf{H}}}_{\mathrm{cp}}\right),$$

(15)

where

$${\widehat{\mathbf{H}}}_{\mathrm{c}}=\sum _{\alpha }|\alpha \rangle {ϵ}_{\alpha }\langle \alpha |$$

(16)

is the one-electron term,

$${\widehat{\mathbf{H}}}_{\mathrm{p}}=\sum _{\mathbf{q}}{ϵ}_{\mathbf{q}}{\widehat{b}}_{\mathbf{q}}^{†}{\widehat{b}}_{\mathbf{q}}$$

(17)

describes the system of noninteracting phonons with wavevector $\mathbf{q}$ and energy ${ϵ}_{\mathbf{q}}$, and

$${\widehat{\mathbf{H}}}_{\mathrm{cp}}=\sum _{\alpha {\alpha }^{\prime }}|\alpha \rangle \left[\sum _{\mathbf{q}}\left(g_{\alpha {\alpha }^{\prime }}^{\mathbf{q},-}{\widehat{b}}_{\mathbf{q}}+g_{{\alpha }^{\prime }\alpha }^{\mathbf{q},+}{\widehat{b}}_{\mathbf{q}}^{†}\right)\right]\langle {\alpha }^{\prime }|$$

(18)

is the (one-electron) carrier–phonon interaction term. Here, ${ϵ}_{\alpha }$ is the energy level corresponding to the generic one-electron eigenstate $|\alpha \rangle$; − and + refer, respectively, to phonon absorption and emission; while the explicit form of the coupling matrix elements $g_{\alpha {\alpha }^{\prime }}^{\mathbf{q},\pm }=g_{\alpha {\alpha }^{\prime }}^{\mathbf{q},\mp *}$ depends on the particular phonon branch (acoustic, optical, etc.) as well as on the coupling mechanism considered (deformation potential, polar coupling, etc.).

Adopting the one-electron picture in Equation (14) and employing the explicit form of the global (carrier plus phonon) Hamiltonian in Equation (15), it is possible to derive the following scattering superoperator acting on the one-electron density matrix ${\rho }_{{\alpha }_1{\alpha }_2}$ [7]:

$${\left.\frac{d{\rho }_{{\alpha }_1{\alpha }_2}}{dt}\right|}_{\mathrm{inco}}=\frac{1}{2}\sum _{{\alpha }_1^{\prime }{\alpha }_2^{\prime }}\left({\mathcal{P}}_{{\alpha }_1{\alpha }_2,{\alpha }_1^{\prime }{\alpha }_2^{\prime }}^{cp}{\rho }_{{\alpha }_1^{\prime }{\alpha }_2^{\prime }}-{\mathcal{P}}_{{\alpha }_1^{\prime }{\alpha }_1^{\prime },{\alpha }_1{\alpha }_2^{\prime }}^{cp*}{\rho }_{{\alpha }_2^{\prime }{\alpha }_2}\right)+\mathrm{H}.\mathrm{c}.$$

(19)

with generalized carrier–phonon scattering rates

$${\mathcal{P}}_{{\alpha }_1{\alpha }_2,{\alpha }_1^{\prime }{\alpha }_2^{\prime }}^{cp}=\sum _{\mathbf{q},\pm }{a}_{{\alpha }_1{\alpha }_1^{\prime }}^{\mathbf{q},\pm }{b}_{{\alpha }_2{\alpha }_2^{\prime }}^{\mathbf{q},\pm *},$$

(20)

where

$$a_{\alpha {\alpha }^{\prime }}^{\mathbf{q},\pm }=\sqrt{\frac{2\pi \left(n_{\mathbf{q}}^{°}+\frac{1}{2}\pm \frac{1}{2}\right)}{ħ}}{g}_{\alpha {\alpha }^{\prime }}^{\mathbf{q},\pm },b_{\alpha {\alpha }^{\prime }}^{\mathbf{q},\pm }=a_{\alpha {\alpha }^{\prime }}^{\mathbf{q},\pm }\delta \left({ϵ}_{\alpha }-{ϵ}_{{\alpha }^{\prime }}\pm {ϵ}_{\mathbf{q}}\right),$$

(21)

and $n_{\mathbf{q}}^{°}$ denotes the equilibrium phonon distribution.

It is worth stressing that the diagonal elements (${\alpha }_1{\alpha }_1^{\prime }={\alpha }_2{\alpha }_2^{\prime }$) of the generalized rates in Equation (20) coincide once again with the usual Fermi’s-golden-rule prescription applied to carrier–phonon interaction [7]:

$$P_{\alpha \to {\alpha }^{\prime }}\equiv {\mathcal{P}}_{{\alpha }^{\prime }{\alpha }^{\prime },\alpha \alpha }^{cp}=\frac{2\pi }{ħ}\sum _{\mathbf{q},\pm }{\left|g_{{\alpha }^{\prime }\alpha }^{\mathbf{q},\pm }\right|}^2\delta \left({ϵ}_{{\alpha }^{\prime }}-{ϵ}_{\alpha }\pm {ϵ}_{\mathbf{q}}\right).$$

(22)

To investigate the impact of possible positivity violations induced by the conventional (i.e., non-Lindblad) adiabatic treatment reviewed so far, we consider a prototypical quantum-dot nanostructure described again via the simple two-level system introduced in Section 1. The latter is characterized by two electronic states only ($\alpha \equiv \{a,b\}$) with an energy splitting ${\Delta }_{\mathrm{c}}={ϵ}_b-{ϵ}_a$ and is described by the two-by-two density matrix in Equation (1). Moreover, regardless of the specific phonon mode considered, we shall adopt $\mathbf{q}$-independent coupling matrices of the following form:

$$\left(\begin{array}{cc}{g}_{bb}^{\mathbf{q},\pm }& g_{ba}^{\mathbf{q},\pm }\\ g_{ab}^{\mathbf{q},\pm }& g_{aa}^{\mathbf{q},\pm }\end{array}\right)=\left(\begin{array}{cc}0& g\\ g& 0\end{array}\right).$$

(23)

This corresponds to neglecting diagonal coupling terms ($a\to a$ and $b\to b$); the related phonon-induced energy renormalizations are relatively small in conventional solid-state systems and have negligible impact on energy dissipation and decoherence phenomena.

As anticipated in Section 1, in this particular case, the basic positivity requirement amounts to asking that whether the eigenvalues ${\mathsf{\Lambda }}_{\pm }$ in Equation (4) or, equivalently, the determinant $f_a{f}_b-{\left|p\right|}^2$ are nonnegative. In order to check the positivity of the density matrix, its eigenvalue analysis is therefore imperative; more specifically, since ${\mathsf{\Lambda }}_{+}\ge {\mathsf{\Lambda }}_{-}$, it is enough to check the nonnegativity of the eigenvalue ${\mathsf{\Lambda }}_{-}$.

For the simulated experiments presented here below, we are assuming as an initial condition a low-density Bell state, namely $f_a\left(0\right)=f_b\left(0\right)=p\left(0\right)\ll 1$, and investigating energy dissipation and decoherence induced on the two-level system by an acoustic-like phonon mode, i.e., a linear-dispersion mode characterized by a bandwidth $0\le {ϵ}_{\mathbf{q}}\le {\Delta }_{\mathrm{p}}$ much greater than the interlevel splitting ${\Delta }_{\mathrm{c}}$.

Figure 2 shows the results obtained via the Markovian dissipation model in Equation (19) equipped with the non-Lindblad scattering rates in Equation (20) in the low-density and low-temperature limit for ${\Delta }_{\mathrm{c}}=4$ meV and for a relatively strong carrier–phonon interaction. More specifically, in order to mimic carrier-acoustic phonon scattering in GaN-based nanomaterials, the coupling coefficient g in Equation (23) and the phonon velocity have been chosen such to produce a semiclassical scattering time $P_{b\to a}^{-1}=0.3$ ps (see Equation (22)), which corresponds to an effective interlevel coupling energy ${\Delta }_{\mathrm{cp}}\equiv ħP_{b\to a}$ of about 2 meV. For the above system parameters, we deal with a dimensionless coupling coefficient $\eta ={\Delta }_{\mathrm{cp}}/{\Delta }_{\mathrm{c}}=0.5$. Such a relatively strong carrier–phonon coupling is indeed typical of GaN-based quantum-dot structures. As we can see, in spite of the typical energy relaxation versus decoherence scenario reported in the upper (population) and middle (polarization) panels, here, the lower panel unequivocally displays negative values (though extremely small) of the density-matrix eigenvalue ${\mathsf{\Lambda }}_{-}$.

To show that the small/negligible positivity violations reported in Figure 2 may also become much larger and thus highly problematic, we have repeated the simulated experiment in Figure 2, reducing the interlevel energy splitting by a factor of four, ${\Delta }_{\mathrm{c}}=1$ meV, which amounts to increasing the coupling coefficient $\eta$ up to 2. The results are reported in Figure 3. As one can see, in spite of a very similar population dynamics (compare upper panels in Figure 2 and Figure 3), we deal with a significant slowdown of the polarization decay (middle panels). Moreover the density-matrix eigenvalue ${\mathsf{\Lambda }}_{-}$ is negative (lower panel) and is an absolute value, now nearly comparable with the typical population and polarization ones. The new simulated experiment in Figure 3 constitutes a non-ambiguous proof that positivity violations, typical of non-Lindblad Markov models, may constitute a severe limitation in the study of solid-state quantum materials and related devices characterized by a strong carrier–phonon coupling, like, e.g., GaN-based quantum-dot nanostructures.

A closer inspection shows that the physical origin of the above positivity violations is an anomalous underestimation of decoherence processes, leading to a slowdown of the polarization decay. This, in turn, leads to negative values of the density-matrix determinant ($f_a{f}_b-{\left|p\right|}^2$) and thus to negative eigenvalues. This limitation is indeed a peculiar feature of the non-Lindblad Markov treatments reviewed so far, and can be easily avoided adopting the alternative (Lindblad-type) adiabatic-decoupling scheme presented below.

2.2. Lindblad-Type Adiabatic-Decoupling Scheme

As anticipated, the choice of the adiabatic-decoupling scheme is definitely not unique. Compared to the conventional approach recalled so far, the alternative adiabatic scheme proposed in Reference [19], based on a time symmetrization between microscopic and macroscopic scales, enables one to replace the non-Lindblad incoherent contribution in Equation (6) with the following Lindblad superoperator [7]:

$${\left.\frac{d\widehat{\rho }}{dt}\right|}_{\mathrm{inco}}=\sum _s\left({\widehat{\mathbf{A}}}^s\widehat{\rho }{\widehat{\mathbf{A}}}^{s†}-\frac{1}{2}\left\{{\widehat{\mathbf{A}}}^{s†}{\widehat{\mathbf{A}}}^s,\widehat{\rho }\right\}\right)$$

(24)

with

$${\widehat{\mathbf{A}}}^s=\underset{\overline{ϵ}\to 0}{lim}{\left(\frac{2{\overline{ϵ}}^2}{\pi {ħ}^6}\right)}^{\frac{1}{4}}{\int }_{-\infty }^{+\infty }{e}^{-\frac{{\widehat{\mathbf{H}}}_{°}{t}^{\prime }}{iħ}}{\widehat{\mathbf{H}}}_s^{\prime }{e}^{\frac{{\widehat{\mathbf{H}}}_{°}{t}^{\prime }}{iħ}}{e}^{-{\left(\frac{\overline{ϵ}{t}^{\prime }}{ħ}\right)}^2}d{t}^{\prime }.$$

(25)

As discussed in Reference [19], here, $\overline{ϵ}$ plays the role of an energy broadening induced by a finite collision duration and/or to a finite single-particle life-time [19].

The new dissipation model in Equation (24) can still be expressed via the global density-matrix superoperator in Equation (8) provided to replace the non-Lindblad rates in Equation (9) with the following (Lindblad-type) version:

$${\mathcal{P}}_{i_1{i}_2,i_1^{\prime }{i}_2^{\prime }}^s={\mathbf{A}}_{i_1{i}_1^{\prime }}^s{\mathbf{A}}_{i_2{i}_2^{\prime }}^{s*}.$$

(26)

While their diagonal (i.e., semiclassical) elements ($i_1{i}_1^{\prime }=i_2{i}_2^{\prime }$) coincide again with the standard Fermi’s-golden-rule prescription in Equation (10), the new (Lindblad) version in Equation (26) exhibits a more symmetric structure, a clear fingerprint of the time symmetrization previously mentioned.

Starting from the Lindblad-type global superoperator in Equation (24) and adopting once again the carrier–phonon mean-field approximation via the factorization scheme in Equation (12), it is possible to get an effective dissipation model for the reduced density-matrix operator ${\widehat{\rho }}_{\mathrm{c}}$, which is still of Lindblad type [7]:

$${\left.\frac{d{\widehat{\rho }}_{\mathrm{c}}}{dt}\right|}_{\mathrm{inco}}=\sum _s\left({\widehat{\mathbf{A}}}_{\mathrm{c}}^s{\widehat{\rho }}_{\mathrm{c}}{\widehat{\mathbf{A}}}_{\mathrm{c}}^{s†}-\frac{1}{2}\left\{{\widehat{\mathbf{A}}}_{\mathrm{c}}^{s†}{\widehat{\mathbf{A}}}_{\mathrm{c}}^s,{\widehat{\rho }}_{\mathrm{c}}\right\}\right).$$

(27)

The explicit form of the new electronic operators ${\widehat{\mathbf{A}}}_{\mathrm{c}}^s$ can be derived starting from the global Lindblad operators ${\widehat{\mathbf{A}}}^s$ in Equation (25). Once again, in spite of their very same formal structure, Equations (24) and (27) describe the system dynamics at different levels.

Adopting once again the one-electron picture previously introduced (see Equation (14)) and applying the Lindblad treatment recalled so far to the case of the carrier-plus-phonon system in Equation (15), it is possible to get the very same density-matrix superoperator in Equation (19) provided to replace the non-Lindblad carrier–phonon rates in Equation (20) with the following Lindblad-type version [7]:

$${\mathcal{P}}_{{\alpha }_1{\alpha }_2,{\alpha }_1^{\prime }{\alpha }_2^{\prime }}=\sum _{\mathbf{q},\pm }{A}_{{\alpha }_1{\alpha }_1^{\prime }}^{\mathbf{q},\pm }{A}_{{\alpha }_2{\alpha }_2^{\prime }}^{\mathbf{q},\pm *},$$

(28)

where

$$A_{\alpha {\alpha }^{\prime }}^{\mathbf{q},\pm }=\sqrt{\frac{2\pi \left(n_{\mathbf{q}}^{°}+\frac{1}{2}\pm \frac{1}{2}\right)}{ħ}}{g}_{\alpha {\alpha }^{\prime }}^{\mathbf{q},\pm }{D}_{\alpha {\alpha }^{\prime }}^{\mathbf{q},\pm }$$

(29)

and

$$D_{\alpha {\alpha }^{\prime }}^{\mathbf{q},\pm }=\underset{\overline{ϵ}\to 0}{lim}\frac{e^{-{\left(\frac{{ϵ}_{\alpha }-{ϵ}_{{\alpha }^{\prime }}\pm {ϵ}_{\mathbf{q}}}{2\overline{ϵ}}\right)}^2}}{{\left(2\pi {\overline{ϵ}}^2\right)}^{\frac{1}{4}}}.$$

(30)

We stress that the diagonal elements (${\alpha }_1{\alpha }_1^{\prime }={\alpha }_2{\alpha }_2^{\prime }$) of the new (Lindblad-type) carrier–phonon rates in Equation (28) coincide once again with the usual Fermi’s-golden-rule prescription in Equation (22).

To concretely test the quality of the alternative adiabatic-decoupling scheme reviewed so far, we have repeated the previous simulated experiments in Figure 2 and Figure 3, replacing the non-Lindblad carrier–phonon scattering rates in Equation (20) with the Lindblad ones in Equation (28). The new results, reported in Figure 4 and Figure 5, fully confirm the absence of positivity violations (see the nonnegative eigenvalue profiles in the lower panels), as expected for any Lindblad-type scattering model. A closer comparison between the former (non-Lindblad) and the present (Lindblad) simulations shows that, in spite of a very similar population dynamics, the Lindblad-type dissipation model is always characterized by a faster polarization decay, thus preventing in any case the positivity violations of the non-Lindblad one reported in Figure 2 and Figure 3.

2.3. Generalization to the Nonlinear Regime

As anticipated, while in the low-density limit, dissipation versus decoherence phenomena can safely be described via the simple one-electron picture previously introduced, at high carrier concentrations, a genuine many-body treatment is imperative [7]; this is typically accomplished via the so-called single-particle picture, based on a mean-field (or Hartree–Fock) treatment of our many-electron system. Indeed, many of the physical quantities of interest in the study of solid-state quantum devices are described via single-particle electronic operators of the following form:

$${\widehat{\mathbf{G}}}_{\mathrm{c}}=\sum _{{\alpha }_1{\alpha }_2}{G}_{{\alpha }_1{\alpha }_2}{\widehat{c}}_{{\alpha }_1}^{†}{\widehat{c}}_{{\alpha }_2},$$

(31)

where ${\widehat{c}}_{\alpha }^{†}$ and ${\widehat{c}}_{\alpha }$ are the usual creation and destruction operators over the electronic single-particle states $|\alpha \rangle$. Recalling that, for any electronic operator, one has $\langle {\mathbf{G}}_{\mathrm{c}}\rangle =\mathrm{tr}\left\{\widehat{\rho }{\widehat{\mathbf{G}}}_{\mathrm{c}}\right\}=\mathrm{tr}{\left\{{\widehat{\rho }}_{\mathrm{c}}{\widehat{\mathbf{G}}}_{\mathrm{c}}\right\}}_{\mathrm{c}}$, the average value of the single-particle operator in Equation (31) can be written as

$$\langle {\mathbf{G}}_{\mathrm{c}}\rangle =\sum _{{\alpha }_1{\alpha }_2}{\rho }_{{\alpha }_1{\alpha }_2}{G}_{{\alpha }_2{\alpha }_1}$$

(32)

where

$${\rho }_{{\alpha }_1{\alpha }_2}=\mathrm{tr}{\left\{{\widehat{c}}_{{\alpha }_2}^{†}{\widehat{c}}_{{\alpha }_1}{\widehat{\rho }}_{\mathrm{c}}\right\}}_{\mathrm{c}}$$

(33)

is the single-particle density matrix. It is worth stressing that, in the low-density limit, the single- particle density matrix in Equation (33) is equivalent to the one-electron density matrix introduced in Equation (14).

To study the time evolution of single-particle quantities, such as total carrier density, mean kinetic energy, charge current, and so on, it is then vital to derive a closed equation of motion for the above single-particle density matrix. Combining its definition in Equation (33) with the many-electron Lindblad dynamics in Equation (27) and employing the cyclic property of the trace, one gets

$${\left.\frac{d{\rho }_{{\alpha }_1{\alpha }_2}}{dt}\right|}_{\mathrm{inco}}=\frac{1}{2}\sum _s\mathrm{tr}{\left\{\left[{\widehat{\mathbf{A}}}_{\mathrm{c}}^{s†},{\widehat{c}}_{{\alpha }_2}^{†}{\widehat{c}}_{{\alpha }_1}\right]{\widehat{\mathbf{A}}}_{\mathrm{c}}^s{\widehat{\rho }}_{\mathrm{c}}\right\}}_{\mathrm{c}}+\mathrm{H}.\mathrm{c}.$$

(34)

To derive a closed equation of motion for the single-particle density matrix, it is now crucial to specify the form/structure of the many-electron Lindblad operators ${\widehat{\mathbf{A}}}_{\mathrm{c}}^s$ in Equation (27), which, in turn, is dictated by the specific interaction mechanism considered.

For the case of a generic carrier–phonon interaction mechanism, the corresponding (one-body) Lindblad operator is always of the following form:

$${\widehat{\mathbf{A}}}_{\mathrm{c}}^s=\sum _{\alpha {\alpha }^{\prime }}{A}_{\alpha {\alpha }^{\prime }}^{\mathrm{cp}}{\widehat{c}}_{\alpha }^{†}{\widehat{c}}_{{\alpha }^{\prime }}.$$

(35)

Equation (35) describes the phonon-induced carrier transition from the initial state ${\alpha }^{\prime }$ to the final state $\alpha$. In this case, the label $s=\mathbf{q},\pm$ corresponds to the emission (+) or absorption (−) of a phonon with wavevector $\mathbf{q}$. By inserting the carrier–phonon Lindblad operator in Equation (35) into Equation (34) and by employing the fermionic anticommutation relations, it is easy to show [21] that the contribution to the system dynamics due to the generic carrier–phonon interaction mechanism ($s=\mathrm{cp}$) involves average values of four fermionic operators of the following form:

$$h_{{\alpha }_3{\alpha }_4,{\alpha }_3^{\prime }{\alpha }_4^{\prime }}=\mathrm{tr}{\left\{{\widehat{c}}_{{\alpha }_3}^{†}{\widehat{c}}_{{\alpha }_4}{\widehat{c}}_{{\alpha }_3^{\prime }}^{†}{\widehat{c}}_{{\alpha }_4^{\prime }}{\widehat{\rho }}_{\mathrm{c}}\right\}}_{\mathrm{c}}.$$

(36)

For the carrier–carrier interaction ($s=\mathrm{cc}$), the Lindblad operator has the general two-body structure:

$${\widehat{\mathbf{A}}}_{\mathrm{c}}^s=\frac{1}{2}\sum _{\alpha \overline{\alpha },{\alpha }^{\prime }{\overline{\alpha }}^{\prime }}{A}_{\alpha \overline{\alpha },{\alpha }^{\prime }{\overline{\alpha }}^{\prime }}^{\mathrm{cc}}{\widehat{c}}_{\alpha }^{†}{\widehat{c}}_{\overline{\alpha }}^{†}{\widehat{c}}_{{\overline{\alpha }}^{\prime }}{\widehat{c}}_{{\alpha }^{\prime }},$$

(37)

The latter describes the transition of the electronic pair from the initial (two-body) state ${\alpha }^{\prime }{\overline{\alpha }}^{\prime }$ to the final state $\alpha \overline{\alpha }$. As shown in Reference [21], by inserting Equation (37) into Equation (34), the contribution to the system dynamics due to carrier–carrier interaction ($s=\mathrm{cc}$) involves average values of eight fermionic operators of the following form:

$$k_{{\alpha }_5{\alpha }_6{\alpha }_7{\alpha }_8,{\alpha }_5^{\prime }{\alpha }_6^{\prime }{\alpha }_7^{\prime }{\alpha }_8^{\prime }}=\mathrm{tr}{\left\{{\widehat{c}}_{{\alpha }_5}^{†}{\widehat{c}}_{{\alpha }_6}^{†}{\widehat{c}}_{{\alpha }_7}{\widehat{c}}_{{\alpha }_8}{\widehat{c}}_{{\alpha }_5^{\prime }}^{†}{\widehat{c}}_{{\alpha }_6^{\prime }}^{†}{\widehat{c}}_{{\alpha }_7^{\prime }}{\widehat{c}}_{{\alpha }_8^{\prime }}{\widehat{\rho }}_{\mathrm{c}}\right\}}_{\mathrm{c}}.$$

(38)

As anticipated, the key step in getting a closed equation of motion for the single-particle density matrix is to adopt the well-known mean-field (or correlation-expansion) approximation [65]; as discussed in Reference [21], employing this approximation scheme and omitting renormalization terms [65], for both carrier–phonon and carrier–carrier scattering, the resulting single-particle equation is given by

$${\left.\frac{d{\rho }_{{\alpha }_1{\alpha }_2}}{dt}\right|}_{\mathrm{inco}}=\frac{1}{2}\sum _{{\alpha }^{\prime }{\alpha }_1^{\prime }{\alpha }_2^{\prime }}\left(\left({\delta }_{{\alpha }_1{\alpha }^{\prime }}-{\rho }_{{\alpha }_1{\alpha }^{\prime }}\right){\mathcal{P}}_{{\alpha }^{\prime }{\alpha }_2,{\alpha }_1^{\prime }{\alpha }_2^{\prime }}^s{\rho }_{{\alpha }_1^{\prime }{\alpha }_2^{\prime }}-\left({\delta }_{{\alpha }^{\prime }{\alpha }_1^{\prime }}-{\rho }_{{\alpha }^{\prime }{\alpha }_1^{\prime }}\right){\mathcal{P}}_{{\alpha }^{\prime }{\alpha }_1^{\prime },{\alpha }_1{\alpha }_2^{\prime }}^{s*}{\rho }_{{\alpha }_2^{\prime }{\alpha }_2}\right)+\mathrm{H}.\mathrm{c}.$$

(39)

with generalized carrier–phonon scattering rates

$${\mathcal{P}}_{{\alpha }_1{\alpha }_2,{\alpha }_1^{\prime }{\alpha }_2^{\prime }}^{s=\mathrm{cp}}=A_{{\alpha }_1{\alpha }_1^{\prime }}^{\mathrm{cp}}{A}_{{\alpha }_2{\alpha }_2^{\prime }}^{\mathrm{cp}*}$$

(40)

and generalized carrier–carrier scattering rates

$${\mathcal{P}}_{{\alpha }_1{\alpha }_2,{\alpha }_1^{\prime }{\alpha }_2^{\prime }}^{s=\mathrm{cc}}=2\sum _{{\overline{\alpha }}_1{\overline{\alpha }}_2,{\overline{\alpha }}_1^{\prime }{\overline{\alpha }}_2^{\prime }}\left({\delta }_{{\overline{\alpha }}_2{\overline{\alpha }}_1}-{\rho }_{{\overline{\alpha }}_2{\overline{\alpha }}_1}\right){\mathcal{A}}_{{\alpha }_1{\overline{\alpha }}_1,{\alpha }_1^{\prime }{\overline{\alpha }}_1^{\prime }}^{\mathrm{cc}}{\mathcal{A}}_{{\alpha }_2{\overline{\alpha }}_2,{\alpha }_2^{\prime }{\overline{\alpha }}_2^{\prime }}^{\mathrm{cc}*}{\rho }_{{\overline{\alpha }}_1^{\prime }{\overline{\alpha }}_2^{\prime }},$$

(41)

where

$${\mathcal{A}}_{\alpha \overline{\alpha },{\alpha }^{\prime }{\overline{\alpha }}^{\prime }}^{\mathrm{cc}}=\frac{1}{4}\left(A_{\alpha \overline{\alpha },{\alpha }^{\prime }{\overline{\alpha }}^{\prime }}^{\mathrm{cc}}-A_{\overline{\alpha }\alpha ,{\alpha }^{\prime }{\overline{\alpha }}^{\prime }}^{\mathrm{cc}}-A_{\alpha \overline{\alpha },{\overline{\alpha }}^{\prime }{\alpha }^{\prime }}^{\mathrm{cc}}+A_{\overline{\alpha }\alpha ,{\overline{\alpha }}^{\prime }{\alpha }^{\prime }}^{\mathrm{cc}}\right)$$

(42)

denote the totally antisymmetric parts of the two-body coefficients in Equation (37).

We stress that, opposite to the carrier–phonon rates in Equation (40), the generalized carrier–carrier rates in (41) are themselves a function of the single-particle density matrix; this is a clear fingerprint of the two-body nature of the carrier–carrier interaction (see below).

The single-particle scattering superoperator in Equation (39) is the result of positive-like (in-scattering) and negative-like (out-scattering) contributions, which are nonlinear functions of the single-particle density matrix. Indeed, in the semiclassical limit [7], namely ${\rho }_{{\alpha }_1{\alpha }_2}=f_{{\alpha }_1}{\delta }_{{\alpha }_1{\alpha }_2}$, the density-matrix equation Equation (39) assumes the expected nonlinear Boltzmann-type form:

$${\left.\frac{d{f}_{\alpha }}{dt}\right|}_{\mathrm{inco}}=\sum _{{\alpha }^{\prime }}\left((1-f_{\alpha })P_{\alpha {\alpha }^{\prime }}^s{f}_{{\alpha }^{\prime }}-(1-f_{{\alpha }^{\prime }})P_{{\alpha }^{\prime }\alpha }^s{f}_{\alpha }\right)$$

(43)

with semiclassical carrier–phonon scattering rates

$$P_{\alpha {\alpha }^{\prime }}^{s=\mathrm{cp}}={\mathcal{P}}_{\alpha \alpha ,{\alpha }^{\prime }{\alpha }^{\prime }}^{s=\mathrm{cp}}={\left|A_{\alpha {\alpha }^{\prime }}^{\mathrm{cp}}\right|}^2$$

(44)

and semiclassical carrier–carrier scattering rates

$$P_{\alpha {\alpha }^{\prime }}^{s=\mathrm{cc}}={\mathcal{P}}_{\alpha \alpha ,{\alpha }^{\prime }{\alpha }^{\prime }}^{s=\mathrm{cc}}=2\sum _{\overline{\alpha }{\overline{\alpha }}^{\prime }}\left(1-f_{\overline{\alpha }}\right){\left|{\mathcal{A}}_{\alpha \overline{\alpha },{\alpha }^{\prime }{\overline{\alpha }}^{\prime }}^{\mathrm{cc}}\right|}^2{f}_{{\overline{\alpha }}^{\prime }}.$$

(45)

The above semiclassical limit clearly shows that the nonlinearity factors $({\delta }_{{\alpha }_1{\alpha }_2}-{\rho }_{{\alpha }_1{\alpha }_2})$ in Equation (39) as well as in Equation (41) can be regarded as the quantum-mechanical generalization of the Pauli factors $(1-f_{\alpha })$ of the conventional Boltzmann theory.

A closer inspection of Equations (39) and (41) as well as of their semiclassical counterparts in Equations (43) and (45) confirms the two-body nature of the carrier–carrier interaction. Indeed, differently from the carrier–phonon scattering, in this case, the density-matrix equation describes the time evolution of a so-called “main carrier” $\alpha$ interacting with a so-called “partner carrier” $\overline{\alpha }$.

Let us finally face the key issue related to the present single-particle treatment, namely the positivity analysis of the nonlinear density-matrix equation in Equation (39). Indeed, for a physical state, the eigenvalues of the single-particle density matrix ${\rho }_{{\alpha }_1{\alpha }_2}$ are necessarily nonnegative but also smaller than one (i.e., Pauli exclusion principle); to preserve such physical nature, it is vital to verify that the scattering-induced time evolution maintains the values of the eigenvalues within the interval [0, 1].

To this end, let us start from the case of carrier–phonon interaction previously discussed, for which the nonlinear equation in Equation (39) (equipped with the generalized rates in Equation (40)) can also be written in a more compact way via the one-electron operators:

$$\widehat{\rho }=\sum _{{\alpha }_1{\alpha }_2}|{\alpha }_1\rangle {\rho }_{{\alpha }_1{\alpha }_2}\langle {\alpha }_2|$$

(46)

and

$$\widehat{A}=\sum _{{\alpha }_1{\alpha }_2}|{\alpha }_1\rangle A_{{\alpha }_1{\alpha }_2}^{\mathrm{cp}}\langle {\alpha }_2|$$

(47)

as

$${\left.\frac{d\widehat{\rho }}{dt}\right|}_{\mathrm{inco}}=\frac{1}{2}\left((\widehat{\mathcal{I}}-\widehat{\rho })\widehat{A}\widehat{\rho }{\widehat{A}}^{†}-{\widehat{A}}^{†}(\widehat{\mathcal{I}}-\widehat{\rho })\widehat{A}\widehat{\rho }\right)+\mathrm{H}.\mathrm{c}.,$$

(48)

where $\widehat{\mathcal{I}}$ denotes the identity operator of the one-electron Hilbert space. We stress that, due to the quantum-mechanical Pauli factors $(\widehat{\mathcal{I}}-\widehat{\rho })$, the above scattering superoperator in Equation (48) is nonlinear in $\widehat{\rho }$ and thus is definitely non-Lindblad. Only in the low-density limit, i.e., $\widehat{\mathcal{I}}-\widehat{\rho }\to \widehat{\mathcal{I}}$, the nonlinear equation in Equation (48) reduces to the Lindblad superoperator

$${\left.\frac{d\widehat{\rho }}{dt}\right|}_{\mathrm{inco}}=\widehat{A}\widehat{\rho }{\widehat{A}}^{†}-\frac{1}{2}\left\{{\widehat{A}}^{†}\widehat{A},\widehat{\rho }\right\},$$

(49)

and the positive-definite character of $\widehat{\rho }$ is thereby ensured.

For high-density conditions, in contrast, no straightforward conclusion can be drawn about the positive-definite character of the corresponding time evolution. Nevertheless, it is easy to show [21] that the nonlinear single-particle equation in Equation (39) does preserve the positive-definite character of $\widehat{\rho }$. To this end, let us consider the instantaneous (i.e., time-dependent) eigenvalues ${\mathsf{\Lambda }}_{\lambda }$ and eigenvectors $|\lambda \rangle$ of the density-matrix operator, i.e.,

$$\widehat{\rho }|\lambda \rangle ={\mathsf{\Lambda }}_{\lambda }|\lambda \rangle ,$$

(50)

which implies that

$${\mathsf{\Lambda }}_{\lambda }=\langle \lambda |\widehat{\rho }|\lambda \rangle .$$

(51)

As anticipated, the eigenvalues in Equation (50) corresponding to a physical state are necessarily nonnegative as well as smaller than one. The scattering-induced time evolution should therefore maintain the values of the ${\mathsf{\Lambda }}_{\lambda }$ within the physical interval [0, 1]; this can be verified by studying the time derivative of the generic eigenvalue in Equation (51):

$$\frac{d{\mathsf{\Lambda }}_{\lambda }}{dt}=\frac{d\langle \lambda |}{dt}\widehat{\rho }|\lambda \rangle +\langle \lambda |\frac{d\widehat{\rho }}{dt}|\lambda \rangle +\langle \lambda |\widehat{\rho }\frac{d|\lambda \rangle }{dt}.$$

(52)

Thanks to the completeness of the basis set $\left\{\right|\lambda \rangle \}$, the time derivative in Equation (52) can also be written as follows:

$$\frac{d{\mathsf{\Lambda }}_{\lambda }}{dt}=\sum _{{\lambda }^{\prime }}\frac{d\langle \lambda |}{dt}|{\lambda }^{\prime }\rangle \langle {\lambda }^{\prime }|\widehat{\rho }|\lambda \rangle +\langle \lambda |\frac{d\widehat{\rho }}{dt}|\lambda \rangle +\sum _{{\lambda }^{\prime }}\langle \lambda |\widehat{\rho }|{\lambda }^{\prime }\rangle \langle {\lambda }^{\prime }|\frac{d|\lambda \rangle }{dt}.$$

Recalling that

$$\langle \lambda |\widehat{\rho }|{\lambda }^{\prime }\rangle ={\mathsf{\Lambda }}_{\lambda }{\delta }_{\lambda {\lambda }^{\prime }},$$

(53)

the result in Equation (53) reduces to

$$\frac{d{\mathsf{\Lambda }}_{\lambda }}{dt}={\mathsf{\Lambda }}_{\lambda }\frac{d\langle \lambda |}{dt}|\lambda \rangle +\langle \lambda |\frac{d\widehat{\rho }}{dt}|\lambda \rangle +{\mathsf{\Lambda }}_{\lambda }\langle \lambda |\frac{d|\lambda \rangle }{dt}.$$

(54)

Taking into account that $\frac{d\langle \lambda |\lambda \rangle }{dt}=0,$ the first and third terms on the right hand side of Equation (54) cancel out exactly, and one finally gets

$$\frac{d{\mathsf{\Lambda }}_{\lambda }}{dt}=\frac{d{\rho }_{\lambda \lambda }}{dt}.$$

(55)

This tells us that the time variation of the eigenvalues ${\mathsf{\Lambda }}_{\lambda }$ coincides with the time variation of the diagonal elements ${\rho }_{\lambda \lambda }$ of the operator $\widehat{\rho }$ within the instantaneous eigenbasis $\left\{\right|\lambda \rangle \}$.

To study the above time derivative, the key step is to examine the explicit form of the single-particle scattering superoperator written in the density-matrix eigenbasis in Equation (50). Taking into account that the generic density-matrix equation in Equation (39) is basis-independent, by replacing the original single-particle basis $\left\{\right|\alpha \rangle \}$ with the density-matrix eigenbasis $\left\{\right|\lambda \rangle \}$ and making use of Equation (53), its diagonal elements turn out to be as follows:

$$\frac{d{\rho }_{\lambda \lambda }}{dt}=\sum _{{\lambda }^{\prime }}\left[(1-{\mathsf{\Lambda }}_{\lambda })P_{\lambda {\lambda }^{\prime }}^s{\mathsf{\Lambda }}_{{\lambda }^{\prime }}-(1-{\mathsf{\Lambda }}_{{\lambda }^{\prime }})P_{{\lambda }^{\prime }\lambda }^s{\mathsf{\Lambda }}_{\lambda }\right],$$

(56)

where

$$P_{\lambda {\lambda }^{\prime }}^s={\mathcal{P}}_{\lambda \lambda ,{\lambda }^{\prime }{\lambda }^{\prime }}^s$$

(57)

are positive-definite quantities given by the diagonal elements of the generalized scattering rates (see Equations (40) and (41)) written in our instantaneous density-matrix eigenbasis. By inserting this last result into Equation (55), we finally get

$$\frac{d{\mathsf{\Lambda }}_{\lambda }}{dt}=\sum _{{\lambda }^{\prime }}\left[(1-{\mathsf{\Lambda }}_{\lambda })P_{\lambda {\lambda }^{\prime }}^s{\mathsf{\Lambda }}_{{\lambda }^{\prime }}-(1-{\mathsf{\Lambda }}_{{\lambda }^{\prime }})P_{{\lambda }^{\prime }\lambda }^s{\mathsf{\Lambda }}_{\lambda }\right].$$

(58)

This result is highly nontrivial: it states that, in spite of the partially coherent nature of the carrier dynamics in Equation (39), the time evolution of the eigenvalues ${\mathsf{\Lambda }}_{\lambda }$ is governed by a nonlinear Boltzmann-type equation, formally identical to the semiclassical result in Equation (43).

We are then finally able to state that the physical interval [0, 1] is the only possible variation range for the eigenvalues ${\mathsf{\Lambda }}_{\lambda }$. Indeed, it is easy to show [21] that, when the latter approach the extremal values, 0 or 1, their time derivatives do not allow them to exit the interval. In particular, a closer inspection of the Boltzmann-like equation in Equation (58) shows the following:

(i)

if one of the eigenvalues ${\mathsf{\Lambda }}_{\lambda }$ is equal to zero, the corresponding time derivative in Equation (58) is always nonnegative;

(ii)

if one of the eigenvalues ${\mathsf{\Lambda }}_{\lambda }$ is equal to one, its time derivative in Equation (58) is always nonpositive.

This leads us to the important conclusion that, for both carrier–phonon and carrier–carrier scattering, the nonlinear single-particle equation in Equation (39) preserves the positive-definite character of the single-particle density matrix.

3. Non-Markovian Dissipation Models

To investigate phonon-induced electronic dissipation versus decoherence via non-Markovian approaches, we shall consider a many-body Hamiltonian of the following form:

$$\widehat{\mathbf{H}}={\widehat{\mathbf{H}}}_{\mathrm{c}}+{\widehat{\mathbf{H}}}_{\mathrm{p}}+{\widehat{\mathbf{H}}}_{\mathrm{cp}},$$

(59)

where

$${\widehat{\mathbf{H}}}_{\mathrm{c}}=\sum _{\alpha }{ϵ}_{\alpha }{\widehat{c}}_{\alpha }^{†}{\widehat{c}}_{\alpha }$$

(60)

is the many-body version of the one-electron Hamiltonian in Equation (16),

$${\widehat{\mathbf{H}}}_{\mathrm{p}}=\sum _{\mathbf{q}}{ϵ}_{\mathbf{q}}{\widehat{b}}_{\mathbf{q}}^{†}{\widehat{b}}_{\mathbf{q}}$$

(61)

is again the phononic Hamiltonian introduced in Section 2.2 (describing the set of noninteracting phonons of wavevector $\mathbf{q}$ and energy ${ϵ}_{\mathbf{q}}$), while

$${\widehat{\mathbf{H}}}_{\mathrm{cp}}=\sum _{\alpha {\alpha }^{\prime },\mathbf{q}}\left(g_{\alpha {\alpha }^{\prime }}^{\mathbf{q},-}{\widehat{c}}_{\alpha }^{†}{\widehat{b}}_{\mathbf{q}}{\widehat{c}}_{{\alpha }^{\prime }}+g_{\alpha {\alpha }^{\prime }}^{\mathbf{q},+}{\widehat{c}}_{{\alpha }^{\prime }}^{†}{\widehat{b}}_{\mathbf{q}}^{†}{\widehat{c}}_{\alpha }\right)$$

(62)

is the many-body version of the one-electron carrier–phonon Hamiltonian in Equation (18).

According to the general quantum-kinetic treatment first proposed in Reference [27], later reviewed in Reference [65], and then extended in Reference [44], we now introduce the set of kinetic variables. In particular, the carrier subsystem is described by the single-particle density matrix (see also Equation (33))

$${\rho }_{{\alpha }_1{\alpha }_2}=\langle {\widehat{c}}_{{\alpha }_2}^{†}{\widehat{c}}_{{\alpha }_1}\rangle$$

(63)

and the phonon subsystem is described by the coherent-phonon amplitude $B_{\mathbf{q}}=\langle {\widehat{b}}_{\mathbf{q}}\rangle$ and by the phonon correlation function

$$n_{\mathbf{q}}=\langle {\widehat{b}}_{\mathbf{q}}^{†}{\widehat{b}}_{\mathbf{q}}\rangle -B_{\mathbf{q}}^{*}{B}_{\mathbf{q}},$$

(64)

where $\langle \dots \rangle =\mathrm{tr}\{\dots \widehat{\rho }\}$.

Starting from the usual Heisenberg equations of motion for the carrier ${\widehat{c}}_{\alpha }$ and phonon ${\widehat{b}}_{\mathbf{q}}$ operators, corresponding to the total many-body Hamiltonian in Equation (59), the dynamics of the coupled carrier–phonon system is governed by the following set of coupled equations:

$$\begin{array}{ccc}\frac{d{\rho }_{{\alpha }_1{\alpha }_2}}{dt}\hfill & =& \frac{1}{ıħ}\left({ϵ}_{{\alpha }_1}-{ϵ}_{{\alpha }_2}\right){\rho }_{{\alpha }_1{\alpha }_2}+\left[\frac{1}{ıħ}\sum _{{\alpha }_3,\mathbf{q}}\left(g_{{\alpha }_1{\alpha }_3}^{\mathbf{q},-}{\varrho }_{{\alpha }_3{\alpha }_2}^{\mathbf{q}}+g_{{\alpha }_3{\alpha }_1}^{\mathbf{q},+}{\varrho }_{{\alpha }_2{\alpha }_3}^{\mathbf{q}*}\right)+\mathrm{H}.\mathrm{c}.\right]\hfill \\ \frac{d{B}_{\mathbf{q}}}{dt}\hfill & =& \frac{1}{ıħ}\left({ϵ}_{\mathbf{q}}{B}_{\mathbf{q}}+\sum _{{\alpha }_1{\alpha }_2}{g}_{{\alpha }_1{\alpha }_2}^{\mathbf{q},+}{\rho }_{{\alpha }_1{\alpha }_2}\right)\hfill \\ \frac{d{n}_{\mathbf{q}}}{dt}\hfill & =& -\left(\frac{1}{ıħ}\sum _{{\alpha }_1{\alpha }_2}{g}_{{\alpha }_1{\alpha }_2}^{\mathbf{q},-}{\varrho }_{{\alpha }_2{\alpha }_1}^{\mathbf{q}}+\mathrm{c}.\mathrm{c}.\right)+\left(\frac{1}{ıħ}\sum _{{\alpha }_1{\alpha }_2}{g}_{{\alpha }_1{\alpha }_2}^{\mathbf{q},-}{\rho }_{{\alpha }_2{\alpha }_1}{B}_{\mathbf{q}}+\mathrm{c}.\mathrm{c}.\right),\hfill \end{array}$$

(65)

where “c.c.” denotes the complex conjugate and

$${\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}=\langle {\widehat{c}}_{{\alpha }_2}^{†}{\widehat{b}}_{\mathbf{q}}{\widehat{c}}_{{\alpha }_1}\rangle ,$$

(66)

are the so-called phonon-assisted density matrices, accounting for the quantum-mechanical phase coherence between carriers and phonons [13,65]. Due to the presence of the latter, the above set of equations is not closed. In particular, the corresponding equation of motion

$$\begin{array}{ccc}\frac{d{\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}}{dt}\hfill & =& \frac{1}{ıħ}\left({ϵ}_{{\alpha }_1}-{ϵ}_{{\alpha }_2}+{ϵ}_{\mathbf{q}}\right){\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}-\frac{1}{ıħ}\sum _{{\alpha }_3,{\mathbf{q}}^{\prime }}\left[g_{{\alpha }_3{\alpha }_2}^{{\mathbf{q}}^{\prime },-}\langle {\widehat{c}}_{{\alpha }_3}^{†}{\widehat{c}}_{{\alpha }_1}{\widehat{b}}_{{\mathbf{q}}^{\prime }}{\widehat{b}}_{\mathbf{q}}\rangle +g_{{\alpha }_2{\alpha }_3}^{{\mathbf{q}}^{\prime },+}\langle {\widehat{c}}_{{\alpha }_3}^{†}{\widehat{c}}_{{\alpha }_1}{\widehat{b}}_{{\mathbf{q}}^{\prime }}^{†}{\widehat{b}}_{\mathbf{q}}\rangle \right.\hfill \\ & -& \left.g_{{\alpha }_1{\alpha }_3}^{{\mathbf{q}}^{\prime },-}\langle {\widehat{c}}_{{\alpha }_2}^{†}{\widehat{c}}_{{\alpha }_3}{\widehat{b}}_{\mathbf{q}}{\widehat{b}}_{{\mathbf{q}}^{\prime }}\rangle -g_{{\alpha }_3{\alpha }_1}^{{\mathbf{q}}^{\prime },+}\langle {\widehat{c}}_{{\alpha }_2}^{†}{\widehat{c}}_{{\alpha }_3}{\widehat{b}}_{\mathbf{q}}{\widehat{b}}_{{\mathbf{q}}^{\prime }}^{†}\rangle \right]+\frac{1}{ıħ}\sum _{{\alpha }_3{\alpha }_4}{g}_{{\alpha }_3{\alpha }_4}^{\mathbf{q},+}\langle {\widehat{c}}_{{\alpha }_2}^{†}{\widehat{c}}_{{\alpha }_4}^{†}{\widehat{c}}_{{\alpha }_3}{\widehat{c}}_{{\alpha }_1}\rangle ,\hfill \end{array}$$

(67)

brings in supplementary, higher-order, kinetic variables given by expectation values of four operators, namely carrier–phonon as well as carrier–carrier two-particle density matrices. The resulting infinite hierarchy of equations may be truncated at some level via a mean-field approximation, assuming that the role played by correlations gets less effective as the number of the involved carriers and/or phonons increases [65]. Within this scheme, the generic phonon-assisted density matrix in Equation (66) is decomposed into its mean-field factorization and a remaining correlation function:

$${\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}={\rho }_{{\alpha }_1{\alpha }_2}{B}_{\mathbf{q}}+\delta {\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}.$$

(68)

At first order, one replaces the phonon-assisted density matrix ${\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}$ in Equation (65) with just its mean-field factorization ${\rho }_{{\alpha }_1{\alpha }_2}{B}_{\mathbf{q}}$. The dynamics of the carrier subsystem then turns out to be described by a fully coherent (coh) term given by

$${\left.\frac{d{\rho }_{{\alpha }_1{\alpha }_2}}{dt}\right|}_{\mathrm{coh}}=\frac{1}{ıħ}\sum _{{\alpha }_3}\left({ϵ}_{{\alpha }_1{\alpha }_3}{\rho }_{{\alpha }_3{\alpha }_2}-{\rho }_{{\alpha }_1{\alpha }_3}{ϵ}_{{\alpha }_3{\alpha }_2}\right),$$

(69)

where

$${ϵ}_{\alpha {\alpha }^{\prime }}={ϵ}_{\alpha }{\delta }_{\alpha {\alpha }^{\prime }}+\left(\sum _{\mathbf{q}}{g}_{\alpha {\alpha }^{\prime }}^{\mathbf{q},-}{B}_{\mathbf{q}}+\mathrm{H}.\mathrm{c}.\right)$$

(70)

indicates the single-particle energy spectrum renormalized by off-diagonal first-order carrier–phonon contributions; in the absence of coherent phonons, $B_{\mathbf{q}}=0$ and such renormalization vanishes. As far as the phonon subsystem is concerned, in the equation of motion for the correlation function $n_{\mathbf{q}}$, the first-order dynamics is absent and

$${\left.\frac{d{n}_{\mathbf{q}}}{dt}\right|}_{\mathrm{coh}}=0,$$

(71)

meaning that $n_{\mathbf{q}}$ remains fixed to its initial value, typically the equilibrium distribution $n_{\mathbf{q}}^{°}$.

To account for incoherent (dissipation versus decoherence) processes, one needs to proceed one step further in the correlation expansion. In particular, the set of kinetic variables has to be enlarged to include also the carrier–phonon correlations $\delta {\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}$ and the related dynamics. The latter is obtained by applying a mean-field approximation to the two-particle (carrier–phonon and carrier–carrier) density-matrices in Equation (67), which are factorized in terms of all possible lowest-order kinetic variables. The resulting set of coupled kinetic equations describing both coherent (coh) and incoherent (inco) carrier–phonon contributions is

$$\begin{array}{ccc}\hfill \frac{d{\rho }_{{\alpha }_1{\alpha }_2}}{dt}& =& {\left.\frac{d{\rho }_{{\alpha }_1{\alpha }_2}}{dt}\right|}_{\mathrm{coh}}+{\left.\frac{d{\rho }_{{\alpha }_1{\alpha }_2}}{dt}\right|}_{\mathrm{inco}}\hfill \\ \hfill \frac{d\delta {\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}}{dt}& =& \sum _{{\alpha }_3{\alpha }_4}{\mathcal{L}}_{{\alpha }_1{\alpha }_2,{\alpha }_3{\alpha }_4}^{\mathbf{q}}\delta {\varrho }_{{\alpha }_3{\alpha }_4}^{\mathbf{q}}+\frac{1}{ıħ}\sum _{{\alpha }_3{\alpha }_4}{g}_{{\alpha }_4{\alpha }_3}^{\mathbf{q},+}\left(n_{\mathbf{q}}+1\right){\rho }_{{\alpha }_4{\alpha }_2}\left({\delta }_{{\alpha }_1{\alpha }_3}-{\rho }_{{\alpha }_1{\alpha }_3}\right)\hfill \\ & -& \frac{1}{ıħ}\sum _{{\alpha }_3{\alpha }_4}{g}_{{\alpha }_4{\alpha }_3}^{\mathbf{q},+}{n}_{\mathbf{q}}{\rho }_{{\alpha }_1{\alpha }_3}\left({\delta }_{{\alpha }_4{\alpha }_2}-{\rho }_{{\alpha }_4{\alpha }_2}\right)\hfill \\ \hfill \frac{d{B}_{\mathbf{q}}}{dt}& =& \frac{1}{ıħ}\left({ϵ}_{\mathbf{q}}{B}_{\mathbf{q}}+\sum _{{\alpha }_1{\alpha }_2}{g}_{{\alpha }_1{\alpha }_2}^{\mathbf{q},+}{\rho }_{{\alpha }_1{\alpha }_2}\right)\hfill \\ \hfill \frac{d{n}_{\mathbf{q}}}{dt}& =& {\left.\frac{d{n}_{\mathbf{q}}}{dt}\right|}_{\mathrm{inco}}\hfill \end{array}$$

(72)

where

$${\left.\frac{d{\rho }_{{\alpha }_1{\alpha }_2}}{dt}\right|}_{\mathrm{inco}}=\frac{1}{ıħ}\sum _{{\alpha }_3,\mathbf{q}}\left(g_{{\alpha }_1{\alpha }_3}^{\mathbf{q},-}\delta {\varrho }_{{\alpha }_3{\alpha }_2}^{\mathbf{q}}+g_{{\alpha }_3{\alpha }_1}^{\mathbf{q},+}\delta {\varrho }_{{\alpha }_2{\alpha }_3}^{\mathbf{q}*}\right)+\mathrm{H}.\mathrm{c}.,$$

(73)

$${\left.\frac{d{n}_{\mathbf{q}}}{dt}\right|}_{\mathrm{inco}}=-\frac{1}{ıħ}\sum _{{\alpha }_1{\alpha }_2}{g}_{{\alpha }_1{\alpha }_2}^{\mathbf{q},-}\delta {\varrho }_{{\alpha }_2{\alpha }_1}^{\mathbf{q}}+\mathrm{c}.\mathrm{c}.,$$

(74)

while

$${\mathcal{L}}_{{\alpha }_1{\alpha }_2,{\alpha }_3{\alpha }_4}^{\mathbf{q}}=\frac{1}{ıħ}\left({ϵ}_{{\alpha }_1{\alpha }_3}{\delta }_{{\alpha }_2{\alpha }_4}-{\delta }_{{\alpha }_1{\alpha }_3}{ϵ}_{{\alpha }_4{\alpha }_2}+{ϵ}_{\mathbf{q}}{\delta }_{{\alpha }_1{\alpha }_3}{\delta }_{{\alpha }_2{\alpha }_4}\right)$$

(75)

is the effective Liouville superoperator acting on the correlation function $\delta {\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}$.

The dynamics of the carrier subsystem, described by the set in Equation (72) is non-Markovian, that is, nonlocal in time. Indeed, by formally integrating the equation of motion for the phonon-assisted correlation function, $\delta {\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}$, and inserting the result into the equation of motion for the electronic density matrix, ${\rho }_{{\alpha }_1{\alpha }_2}$, one gets an integrodifferential equation showing that the time-derivative of the latter at a given time depends on its values at all previous times [65].

The third equation in Equation (72) describes the coherent-phonon dynamics and explicitly shows that, in general, the carrier–phonon coupling itself may generate a finite coherent-phonon amplitude $B_{\mathbf{q}}$ even when the latter is initially null. The excitation of coherent phonons in semiconductors, first observed in the surface field of n-doped GaAs by optical excitation [66], is nowadays a well-established phenomenon and can be realized in systems with sufficiently low symmetry via a properly tailored ultrafast optical excitation.

Finally, the fourth equation in Equation (72), containing the incoherent term in Equation (74), is often referred to as the hot-phonon equation, since it describes how the phonon distribution $n_{\mathbf{q}}$ may be driven out of its thermal-equilibrium value $n_{\mathbf{q}}^{°}$ via the electron–phonon coupling. Such hot-phonon effects are known to play a key role in specifically designed and optimized quantum-cascade devices [67,68,69].

Applying to the coupled set of kinetic equations in Equation (72) of the alternative Markov treatment recalled in Section 2.2 and Section 2.3, the phonon-assisted correlation $\delta {\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}$ can be adiabatically eliminated, allowing for the derivation of the (positive-definite) nonlinear density-matrix equation in Equation (39) equipped with the carrier–phonon scattering rates in Equation (28).

In the low-density limit, Pauli factors can be neglected (${\delta }_{\alpha {\alpha }^{\prime }}-{\rho }_{\alpha {\alpha }^{\prime }}\to {\delta }_{\alpha {\alpha }^{\prime }}$) and hot-phonon effects are absent ($n_{\mathbf{q}}\to n_{\mathbf{q}}^{°}$); the quantum-kinetic set in Equation (72) then assumes a much simpler form. In addition to this, an initially null phonon amplitude ($B_{\mathbf{q}}\left(0\right)=0$) will remain zero at any later time, implying that all first-order carrier–phonon energy renormalizations in Equation (70) vanish (${ϵ}_{\alpha {\alpha }^{\prime }}\to {ϵ}_{\alpha }{\delta }_{\alpha {\alpha }^{\prime }}$) and that the electron–phonon correlation function coincides with the phonon-assisted density matrix ($\delta {\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}\to {\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}$). Under these simplifying low-density conditions, the full set of quantum-kinetic equations in Equation (72) reduces to

$$\begin{array}{ccc}\frac{d{\rho }_{{\alpha }_1{\alpha }_2}}{dt}\hfill & =& \frac{1}{ıħ}\left({ϵ}_{{\alpha }_1}-{ϵ}_{{\alpha }_2}\right){\rho }_{{\alpha }_1{\alpha }_2}+\left[\frac{1}{ıħ}\sum _{{\alpha }_3,\mathbf{q}}\left(g_{{\alpha }_1{\alpha }_3}^{\mathbf{q},-}{\varrho }_{{\alpha }_3{\alpha }_2}^{\mathbf{q}}+g_{{\alpha }_3{\alpha }_1}^{\mathbf{q},+}{\varrho }_{{\alpha }_2{\alpha }_3}^{\mathbf{q}*}\right)+\mathrm{H}.\mathrm{c}.\right]\hfill \\ \frac{d{\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}}{dt}\hfill & =& \frac{1}{ıħ}\left({ϵ}_{{\alpha }_1}-{ϵ}_{{\alpha }_2}+{ϵ}_{\mathbf{q}}\right){\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}+\frac{1}{ıħ}\sum _{{\alpha }_3}\left[g_{{\alpha }_3{\alpha }_1}^{\mathbf{q},+}\left(n_{\mathbf{q}}^{°}+1\right){\rho }_{{\alpha }_3{\alpha }_2}-g_{{\alpha }_2{\alpha }_3}^{\mathbf{q},+}{n}_{\mathbf{q}}^{°}{\rho }_{{\alpha }_1{\alpha }_3}\right].\hfill \end{array}$$

(76)

Applying again the alternative adiabatic-decoupling scheme in Section 2.2 to the low-density coupled set in Equation (76), the latter reduces to the linear scattering superoperator in Equation (19) equipped again with the Lindblad-type carrier–phonon rates in Equation (28).

3.1. Low-Density Analysis

To focus on the comparison between the above recalled quantum-kinetic treatments and their Markovian counterparts, let us consider again the two-level system previously introduced (see Equations (1) and (23)). More specifically, for the simulated experiments presented in this section, the initial condition of the electronic system is still a low-density Bell state, namely $f_a\left(0\right)=f_b\left(0\right)=p\left(0\right)\ll 1$, with null coherent-phonon amplitudes, $B_{\mathbf{q}}\left(0\right)=0$, and phonon-assisted density matrices, ${\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}\left(0\right)=0$.

We start analysing energy dissipation and decoherence effects due to the very same acoustic-like phonon mode considered in Section 2. Figure 6 compares the results of the Lindblad-type Markovian approach (MA) and the quantum-kinetic approach (QKA) for the same cases considered in Figure 2 and Figure 4. Due to the relatively strong carrier–phonon coupling, the quantum-kinetic dynamics shows significant differences from the Markovian one. In particular, the initial time derivative of the QKA excited-level population (upper panel) is null, a well-known feature of non-Markovian models, and the transient dynamics is heavily influenced by energy non-conserving transitions [7]. The combination of these two aspects initially causes a reduction of energy dissipation and decoherence, which however tends to vanish at longer times. The results in the upper and middle panels of Figure 6 undoubtedly show the effectiveness of the QKA in describing a zero-dimensional electronic system strongly coupled to acoustic phonon modes; a typical example, in this respect, is the phase-coherence versus dissipation dynamics in quantum-dot-based nanodevices [3,7]. However, in spite of these facts, positivity violations—albeit rather small—show up when looking at the eigenvalue behaviour (lower panel).

We now move to the analysis of energy dissipation and decoherence effects when optical-like phonon modes are considered. Indeed, in the case of a purely zero-dimensional electronic system coupled to a dispersionless phonon mode (that is, a discrete electron-plus-phonon energy spectrum), the Markov limit is inapplicable and more refined treatments, e.g., based on the polaronic picture [70,71], are necessary. For this reason, we here assume a finite phonon bandwidth ${\Delta }_{\mathrm{p}}$, centered around the electronic interlevel excitation and much smaller than the interlevel splitting ${\Delta }_{\mathrm{c}}$, so that the Markov approach can be still well defined. In particular, in order to mimic electron-optical phonon scattering in GaAs-based nanostructures, we shall employ the following system parameters: ${\Delta }_{\mathrm{c}}=40$ meV, ${\Delta }_{\mathrm{p}}=4$ meV and a semiclassical scattering time $P_{b\to a}^{-1}=0.3$ ps, corresponding to an effective interlevel coupling energy ${\Delta }_{\mathrm{cp}}$ of about 2 meV. The low-density and low-temperature results for this prototypical system are reported in Figure 7, again comparing the Markovian and the non-Markovian schemes. In spite of the weak-coupling regime ($\eta ={\Delta }_{\mathrm{cp}}/{\Delta }_{\mathrm{c}}\simeq 0.05$), the QKA dynamics now differs significantly from the MA counterpart. Indeed, the former displays an almost dissipation-free oscillatory behaviour for the excited-level population and the interlevel polarization, while the latter has the typical exponential decay (upper and middle panels). A more relevant difference, however, is the fact that in the QKA evolution the excited-level population turns negative; this fact matches with the positivity violation in the eigenvalue analysis, reported in the lower panel.

Indeed, in this respect, the relatively small value of the optical-like phonon bandwidth is a critical issue. To investigate this aspect in a deeper way, Figure 8 compares the present QKA population and eigenvalue dynamics (dashed line in the upper and lower panel, respectively), corresponding to ${\Delta }_{\mathrm{p}}=4$ meV, with the results obtained for two different values of the phonon bandwidth. In particular, ${\Delta }_{\mathrm{p}}$ is set to 0 meV (solid line) and 8 meV (dash-dotted line). While, for null ${\Delta }_{\mathrm{p}}$, one gets a cosine-like behavior of the population and the eigenvalue, on increasing ${\Delta }_{\mathrm{p}}$, the negative regions in both the corresponding profiles are significantly reduced and the dissipation effects increased. The quantum-kinetic dynamics displayed in Figure 7 and Figure 8 claim some remarks, mainly on its oscillatory character and positivity violations.

The oscillatory behavior is clear signature of a dissipation-free dynamics; this is a well known effect [28,72] and not an artifact of the QKA: a dispersionless phonon mode, ${\Delta }_{\mathrm{p}}\to 0$, does not induce any dissipation and decoherence in the electronic subsystem. Within the QKA, electronic dissipation versus decoherence phenomena result from a complex interference process involving all different electron and phonon energies. In the weak-coupling limit, all phonon-assisted density-matrix elements are expected to rotate with different frequencies: ${\varrho }_{{\alpha }_1{\alpha }_2}^{\mathbf{q}}\propto exp\left[-ı\left({ϵ}_{{\alpha }_1}-{ϵ}_{{\alpha }_2}+{ϵ}_{\mathbf{q}}\right)/ħ\right]$; this implies that the dynamics of the single-particle density matrix ${\rho }_{{\alpha }_1{\alpha }_2}$ (first equation in Equation (72)) is the outcome of a purely coherent and reversible interference process involving all phonon-assisted correlations. This scenario is analogous to the ultrashort temporal decay of the total polarization in a photoexcited semiconductor [73,74] resulting from the coherent superposition of microscopic polarizations rotating with different frequencies (inhomogeneous broadening), which is present also in the absence of genuine decoherence processes (homogeneous broadening). However, in the case of a zero-dimensional electronic system, the efficacy of the abovementioned interference process requires the presence of a continuum of phonon-assisted density-matrix energies much larger than the typical interstate energy splitting ${ϵ}_{\alpha }-{ϵ}_{{\alpha }^{\prime }}$. This requirement is always fulfilled in the case of the acoustic-phonon mode (see Figure 6). However, the same does not apply to the case of the optical-phonon mode and this is the main origin of the almost dissipation-free results in Figure 7 and Figure 8. It is crucial to stress once again that the anomalous dissipation dynamics pointed out so far is ascribed to the dispersionless-phonon limit only and is by no means related to the weak- versus strong-coupling regime; indeed, all the optical-phonon results presented so far (see Figure 7 and Figure 8) refer to the weak-coupling regime ($\eta \simeq 0.05$).

The presence of unambiguous positivity violation in the QKA results in Figure 7 and Figure 8 represents a severe problem since it arises in the dispersionless regime, ${\Delta }_{\mathrm{p}}\to 0$, where quantum-kinetic approaches are typically invoked as a reliable alternative to—inapplicable—Markov treatments. Therefore, while on the one hand quantum-kinetic approaches have unquestionably proven to explain a variety of ultrafast optical phenomena (like the phonon quantum beats observed, e.g., in Reference [75]) where the key role of uncompleted collisions could not be accounted for within a Markov treatment, on the other hand, the present analysis points out that caution should be used in considering QKA results as faultless overall.

To address this point more extensively, let us start noticing that, in the dispersionless limit (${\Delta }_{\mathrm{p}}\to 0$), each phonon $\mathbf{q}$ is characterized by the very same energy, i.e., ${ϵ}_{\mathbf{q}}\to {ϵ}_b-{ϵ}_a$. The effect of the whole phonon system can be therefore described via a single phonon $\overline{\mathbf{q}}$ resonantly coupled to the two-level system via an effective coupling constant $\overline{g}$. Writing the phonon-assisted density matrix in the following form

$$\left(\begin{array}{cc}{\varrho }_{bb}^{\overline{\mathbf{q}}}& {\varrho }_{ba}^{\overline{\mathbf{q}}}\\ {\varrho }_{ab}^{\overline{\mathbf{q}}}& {\varrho }_{aa}^{\overline{\mathbf{q}}}\end{array}\right)=\left(\begin{array}{cc}{F}_b& P_{ba}\\ P_{ab}& F_a\end{array}\right)$$

(77)

and adopting the usual rotating-wave approximation, the original quantum-kinetic equations in Equation (76) in the low-temperature limit ($n_{\overline{\mathbf{q}}}^{°}=0$) allow us to derive a set of coupled equations for $f_b$ and $P_{ab}$,

$$\begin{array}{ccc}\frac{d{f}_b}{dt}\hfill & =& 2\overline{\omega }\Im \left(P_{ab}\right)\hfill \\ \frac{d{P}_{ab}}{dt}\hfill & =& -ı\overline{\omega }{f}_b,\hfill \end{array}$$

(78)

as well as for p, $F_a$, and $F_b$,

$$\begin{array}{ccc}\frac{dp}{dt}\hfill & =& -ı{\omega }_{°}p-ı\overline{\omega }\left(F_a-F_b\right)\hfill \\ \frac{d{F}_a}{dt}\hfill & =& -ı{\omega }_{°}{F}_a-ı\overline{\omega }p\hfill \\ \frac{d{F}_b}{dt}\hfill & =& -ı{\omega }_{°}{F}_b,\hfill \end{array}$$

(79)

where ${\omega }_{°}={\Delta }_{\mathrm{c}}/ħ$ and $\overline{\omega }=\overline{g}/ħ$. The two sets in Equations (78) and (79) are independent and can be solved analytically; in particular, assuming an initially null phonon-assisted density matrix in Equation (77), the excited-level population evolves in time as

$$f_b\left(t\right)=f_b\left(0\right)cos\left(\sqrt{2}\overline{\omega }t\right).$$

(80)

Such a cosine-like dynamics is in full agreement with the dispersionless-limit behavior depicted in Figure 8, confirming once again the potential limitations of the quantum-kinetic treatment (76) in describing a two-level system coupled to a single phonon $\overline{\mathbf{q}}$.

As a further step, we shall consider the dynamics of the global (electron-plus-one phonon) system introduced above. In particular, the time evolution of its density matrix operator $\widehat{\rho }$ is dictated by the following Liouville–von Neumann equation, corresponding to the total many-body Hamiltonian in Equation (59):

$$\frac{d\widehat{\rho }}{dt}=\frac{1}{ıħ}\left[\widehat{\mathbf{H}},\widehat{\rho }\right].$$

(81)

Figure 9 compares the dynamics of the excited-level population resulting from a numerical solution of the exact global density-matrix equation in Equation (81) and the quantum-kinetic outcome obtained in the case ${\Delta }_{\mathrm{p}}=0$ meV. Both curves show a fully coherent dynamics and the initial zero-derivative behavior, typical of any genuine quantum-mechanical treatment. However, the QKA oscillation period comes out to be a factor of $\sqrt{2}$ larger than the exact-model one and, more importantly, while the exact result is always positive-definite, its quantum-kinetic counterpart is not.

As final step of our low-density analysis, we shall now verify the possible role played by the phonon-bath temperature, which has so far been assumed equal to 0 K. To answer this question, we have therefore repeated the set of simulated experiments in Figure 8 at, e.g., room-temperature. The results are reported in Figure 10 and show a significant suppression of the negative-value regions; nevertheless, also at room temperature, in the dispersionless limit, we deal once again with a dissipation-free cosine-like population dynamics. Its oscillations are faster than the zero-temperature ones and are no more symmetric. Indeed, it is straightforward to show that the finite-temperature generalization of the analytical (zero-temperature) result in Equation (80) is given by

$$f_b\left(t\right)=\tilde{f}+(f_b\left(0\right)-\tilde{f})cos\left(\sqrt{2}\tilde{\omega }t\right)$$

(82)

where

$$\tilde{f}=\left(f_a\left(0\right)+f_b\left(0\right)\right)\frac{{\overline{n}}^{°}}{2{\overline{n}}^{°}+1}$$

(83)

and $\tilde{\omega }=\sqrt{2{\overline{n}}^{°}+1}\overline{\omega }$, where ${\overline{n}}^{°}\equiv n_{\overline{\mathbf{q}}}^{°}$ denotes the thermal occupation number of the dispersionless phonon mode resonantly coupled to our two-level system.

3.2. High-Density Analysis

The goal of this final section is to show if and to which extent the positivity violations encountered in the previous low-density analysis are present at high carrier concentrations as well.

In genuine zero-dimensional systems like, e.g., semiconductor macroatoms [3,7], the presence of a significant quantum-mechanical carrier confinement [76] determines a relevant interlevel energy splitting ${\Delta }_{\mathrm{c}}$ between the two lowest electronic states. In equilibrium or quasiequilibrium conditions, the carrier occupation of all the higher-energy states is therefore negligible, so that the latter may be safely neglected and the electro-optical response of the quantum device may be properly described via our prototypical two-level system. However, opposite to the low-density regime of Section 3.1, in the high-density case, the impact of electronic degeneracy is expected to play an important role. For such a class of quantum systems and related devices, it is then vital to account for Pauli-blocking effects; to this end, the low-density quantum-kinetic set in Equation (76) should be replaced by the original (nonlinear) set in Equation (72).

To investigate the impact of electronic degeneracy, we have repeated the low-temperature simulated experiments of Figure 8, replacing the initial low-density Bell state with the maximally-degenerate Bell state $f_a\left(0\right)=f_b\left(0\right)=p\left(0\right)=\frac{1}{2}$. The results for the excited-level population and eigenvalue dynamics are reported in Figure 11 and clearly show that nonlinear effects have a strong impact on the system time evolution; indeed, comparing the high-density (or maximally-degenerate) results in Figure 11 to the low-density (or non-degenerate) results in Figure 8, we see that the inclusion of the Pauli factors $({\delta }_{\alpha {\alpha }^{\prime }}-{\rho }_{\alpha {\alpha }^{\prime }})$ in Equation (72) leads to a suppression of the positivity-violation signatures in the population dynamics. This is a clear indication that, in such a maximally-degenerate regime, the initial-condition parameter space leading to positivity violations is significantly reduced compared to its low-density counterpart.

4. Summary and Conclusions

The aim of this paper was to provide a cohesive review of the most important dissipation/ decoherence models routinely employed for the design and microscopic simulation of solid-state quantum devices, pointing out reciprocal virtues versus limitations. In particular, we focused on a few critical issues related to conventional Markov models [21] as well as to quantum-kinetic treatments [64], linking them within a common framework. More specifically, thanks to properly designed simulated experiments of a prototypical quantum-dot nanostructure (described via an electronic two-level system coupled to a phonon bath), we reached the following conclusions.

As far as Markov dissipation models are concerned, we have first shown that conventional (i.e., non-Lindblad) adiabatic-decoupling schemes may lead to positivity violations, pointing out regimes where such unphysical behaviours are particularly severe; we have then shown that these limitations may be definitely avoided adopting an alternative (i.e., Lindblad-type) adiabatic-decoupling scheme; the latter has been finally generalized to the nonlinear (i.e., degenerate) regime.

As far as dissipation models based on quantum-kinetic treatments are concerned, we have first shown that, at low temperature and low carrier concentration, the mean-field approximation may lead again to incongruous behaviors, characterized by anomalous decoherence suppression and/or positivity violations; we have then shown that the inclusion of finite-temperature conditions and/or high-density effects leads to a significant reduction of the anomalous behaviors just recalled. It is finally vital to stress that these anomalous behaviours do not affect semiconductor nanostructures with a continuous electron-plus-phonon spectrum; the latter include quantum wells and wires (where the continuous electronic spectrum allows anyway for a proper treatment of dispersionless phonon modes) as well as zero-dimensional electronic systems (quantum dots) coupled to acoustic phonons. In view of the above, the potential limitations of non-Markovian dissipation models studied in this paper (i) are expected to have negligible consequences on previous quantum-kinetic investigations and (ii) may play a significant role in the theoretical modelling of new-generation quantum nanomaterials and nanodevices operating at low carrier density and temperature and characterized by discrete electronic as well as phononic energy spectra. This concerns, in the first instance, conventional semiconductor quantum-dot based nanodevices, particularly in the regime of strong coupling to nearly dispersionless optical modes, as is typically the case for GaN-based nanostructures [77,78].

In summary, the most important conclusion of our investigation is that dissipation-induced positivity violations are ascribed not only to the adiabatic approximation but also to the mean-field approximation; indeed, while in the case of Markov treatments (based on adiabatic as well as mean-field approximations) the problem is easily avoidable by employing Lindblad-type dissipation models, for the case of non-Markovian treatments (based on the mean-field approximation only), a general strategy is still missing.