Ultrafast dynamics of laser-pulse excited semiconductors : non-Markovian quantum kinetic equations with nonequilibrium correlations

Non-Markovian kinetic equations in the second Born approximation are derived for a two-zone semiconductor excited by a short laser pulse. Both collision dynamics and running nonequilibrium correlations are taken into consideration. The energy balance and relaxation of the system to equilibrium are discussed. Results of numerical solution of the kinetic equations for carriers and phonons are presented.


Introduction
Within the last decade, ultrafast irreversible processes induced by laser pulses in semiconductors have attracted considerable attention [1][2][3].Besides a relevant technological interest in microelectronics, the studies of optical and transport properties of semiconductors under high intensity laser radiation provide the motives for the re-examination of some fundamental ideas of quantum kinetics: memory effects, quantum coherence, the effect of many-particle correlations, formation of quasiparticle behavior, etc.It would be beyond the scope of this paper to review the current state of the art in this field (see, e.g., [3]).We will concentrate only on some aspects of the present theories that clarify the aim of our work.The most commonly used thec V.V.Ignatyuk, V.G.Morozov oretical approaches to ultrafast dynamics in semiconductors are the density matrix method (see the recent review in [3]) and the Green's function formalism [2,4,5].In the density matrix method one starts with equations of motion for the single-time carrier distribution functions, polarization, and the phonon distribution function.These equations are not closed; they involve the so-called phonon-assisted density matrices describing correlations between carriers and phonons.On the next level, the equations of motion for the phonon-assisted density matrices contain higher-order correlation functions, so that an infinite hierarchy of coupled equations evolves.The central problem in this approach is the truncation of the hierarchy.Different truncation schemes based on more or less reasonable physical ideas have been proposed [3].It should be noted, however, that the number of equations increases rapidly with the level of truncation.Because of this, it is often necessary to take recourse to further approximations.A weak point in such a strategy is that even quite "reasonable" approximations can lead to non-physical predictions.For instance, some approximations violate the conservation of the total energy for the case of relaxation from a given initial distribution or violate the energy balance for the case of interaction with the exciting light field.Note also that the validity of truncation schemes for extremely ultrashort time scales remains open to question.An alternate approach to quantum kinetics at short time scales is the nonequilibrium Green's function technique [2,6].The basic ingredients of this formalism are two-time quantities: the spectral functions and the correlation functions which satisfy the so-called Kadanoff-Baym equations.In principle, the Kadanoff-Baym equations can serve as a starting point for the derivation of kinetic equations for single-time functions (the Wigner distributions and polarization).To obtain closed kinetic equations, it is necessary to "reconstruct" the two-time correlation functions in terms of their values on the time diagonal, i.e., in terms of the single-time distribution functions.This is a crucial point in the Green's function approach.Formally, the reconstruction problem can be solved exactly by way of very complicated integral equations [7] involving the Wigner functions, the self-energies, and the retarded and advanced Green's functions (or, equivalently, the spectral functions).The dominant terms in these equations, at least for small self-energies, lead to the approximate relation between two-time correlation functions and the Wigner distributions, which is known as the generalized Kadanoff-Baym ansatz (GKBA) [6,7].Non-Markovian kinetic equations derived based on the GKBA are of considerable current use in the theory of ultrafast relaxation processes [2].Unfortunately, the GKBA does not provide a consistent determination of the spectral function.Moreover, spectral functions with non-zero quasiparticle damping lead to non-Markovian kernels in kinetic equations which violate the energy balance at the stage of carrier excitation by the light field and do not conserve total energy at the relaxation stage [2].As a result, the present theory does not yield the correct asymptotic distributions.Some authors proposed spectral functions [6] that at least partly improve the asymptotic behavior of the distribution functions.It should be noted, however, that the attempts to construct an appropriate spectral function in the GKBA for each concrete situation do not provide any insight into the problem of energy conservation in non-Markovian quantum kinetics.A more promising ap-proach is to solve the full two-time Kadanoff-Baym equations where the quasiparticle damping is taken into account self-consistently [6].Because of computer limitations, explicit calculations have so far been restricted to very simple models.In this paper we present a new approach to the theory of ultrafast processes in semiconductors within the framework of non-Markovian quantum kinetics.The basic idea outlined earlier [8,9] is to treat the total energy (or, equivalently, the interaction energy) and the single-particle distributions as independent state parameters.This leads to additional "correlation" terms in non-Markovian collision integrals that ensure the energy balance and yield the correct asymptotic long-time behavior of the distribution functions 1 .The plan of the paper is as follows.In section 2 we outline a formulation of the problem in terms of single-particle distributions.The collision integrals in non-Markovian Born approximation are derived in section 3. We also introduce the time-dependent quasi-temperature which is thermodynamically conjugated to the mean interaction energy and plays a crucial role in our approach.The equation for the quasi-temperature is derived by the method of nonequilibrium statistical operator [11].Section 4 deals with the energy balance equation.It is shown that the total energy is exactly conserved at the relaxation stage, and the explicit expression for the mean value of the nonequilibrium interaction energy is obtained.Relaxation of the carrier and phonon subsystems to thermal equilibrium is discussed in section 5.The problems arising when a nonzero quasiparticle damping is taken into account are discussed in section 6, and results of numerical solution of kinetic equations for carriers and phonons are presented in section 7. The final section contains a discussion and comments on further investigations.

A model of a two-zone laser excited semiconductor
The Hamiltonian of the system will be taken in the form Ĥ(t) = Ĥ0 + Ĥf (t) + Ĥint . (2.1) The kinetic energy operator is given by where a † αk and a αk , α = {c, v}, are creation and annihilation operators for the conduction (α = c) and valence electrons (α = v), and b † q , b q are creation and annihilation operators for phonons.For simplicity, we restrict our consideration to a parabolic two-band semiconductor with where E g is the energy gap, m c is the effective mass of electrons in the conduction band, and m v is the effective mass of holes in the valence band.The field-dependent part of the Hamiltonian, Ĥf (t), describes the interband polarization of the carriers by an external field E(t) and can be written as with the optical matrix elements where u αk is the Bloch function for the band α, and E 0 is the field amplitude vector.Usually, µ αβ (k) are supposed to be constants, µ αβ (k) ≈ d 0 ; they define the Rabi energy ω R = d 0 E 0 of the laser pulse [2,5].The first two terms of the Hamiltonian (2.1) can be presented in the following form: where In real experiments, the time profile of the external electric field can be approximated as a leading edge of the excitation pulse with duration τ pulse and the central frequency ω.For definiteness, we shall take it in the form Since the characteristic wavelength of the light field is about 10 −7 m, which is much larger than the lattice constant, the system may be considered as being spatially homogeneous.The last term, Ĥint , in the Hamiltonian (2.1) describes the carrierphonon interaction and has the standard form with the Fröhlich coupling [2,4] where ω LO is the energy of longitudinal optical (LO) phonons, and * denotes an effective dielectric constant.In what follows, the LO-phonon dispersion is neglected.The carrier density is assumed to be sufficiently low, so that the carrier-carrier Coulomb interaction is of minor importance compared with the interaction with LO-phonons.In this paper, the main emphasis will be done on the derivation of kinetic equations for the distribution functions of carriers, f αβ (k, t), and phonons, n(q, t), which are defined as where we use the notation p(k, t) for the polarization [5].The averages in (2.10) are calculated with a nonequilibrium statistical operator (NSO) ρ(t) which satisfies the von-Neumann equation Calculating the time derivatives of the distribution functions and using the von-Neumann equation, we readily obtain the kinetic equations where [M(k, t), f (k, t)] αβ denotes the matrix elements of the commutator of two matrixes.The collision integrals are given in the general form as (2.13) Equations (2.12) are presented in an exact form with none approximations for the collision integrals.Equations with I αβ (k, t) = 0, I ph (q, t) = 0 correspond to the so-called coherent approximation [2,4].In this approximation, the kinetic equations for the components of the matrix distribution function can be written as These equations are nothing but the Bloch equations which are widely used to describe the optical response of a two-level system.A solution of equations (2.14) for an arbitrary E(t) is not known.Note, however, that these equations yield a useful integral of motion [5] 4|p which shows that the quantities |p(k, t)| 2 and f c (k, t) are generally of the same order.As already discussed, the main source of difficulties in the kinetic theory of ultrafast processes in laser-excited semiconductors is the problem of energy balance at the excitation stage and the energy conservation at the relaxation stage.Following the approach outlined in [8,9], we shall treat the mean interaction energy Ĥint t as an additional quantity characterizing the nonequilibrium state of the system.In this way, the energy balance equation will follow directly from the kinetic equations (2.12) and from the equation of motion for Ĥint t .

Collision integrals and equation for quasi-temperature
To calculate the collision integrals (2.13) explicitly, one has to find a retarded solution of the von-Neumann equation (2.11) as a functional of the state parameters f αβ (k, t ), n(q, t ), and Ĥint t , where t < t.Following the standard procedure (see, e.g., [11]), we first define the auxiliary relevant statistical operator in a generalized Gibbs form, where the partition function is determined by the normalization condition, and the Lagrange multipliers β(t), Λ αβ (k, t), Λ ph (q, t) are determined by the selfconsistency conditions Here and henceforth the symbol • • • t rel stands for averages calculated with ρ rel (t).Clearly the equilibrium statistical operator follows from equation (3.1) as a special case.In thermal equilibrium, β = 1/T eq , where T eq is the temperature, and the Lagrange multipliers conjugated to the distribution functions are given by where µ is the equilibrium chemical potential.Thus, assuming that the evolution of the system starts with thermal equilibrium at time t 0 , one may look for a solution of the von Neumann equation that satisfies the initial condition We shall restrict our consideration to the second-order non-Markovian Born approximation in which the collision integrals (2.13) are explicitly proportional to D 2 (q).Within this approximation, we need a solution of equation (2.11), correct to first order in the interaction H int .This is given by where is the evolution operator with the Hamiltonian (2.6).The symbol exp + {. ..} means the chronologically ordered exponent.Substituting (3.6) into (2.13),we readily obtain the collision integrals in the form The action of U † 0 (t, t ) and U 0 (t, t ) on creation and annihilation operators are calculated explicitly only for phonons: (3.10)For fermions, the corresponding transformations can be written in a symbolic matrix form where we have introduced matrices These time-ordered exponents can be found by solving the matrix equation In general, only numerical solutions are possible.An approximate analytical solution for sufficiently small field amplitudes is considered in Appendix A. Now we are in a position to calculate the first terms in (3.8) and (3.9).We call them the "collision" contributions to the total collision integrals I αβ (k, t), I ph (q, t), in order to distinguish them from the "correlation" contributions associated with the last terms in (3.8) and (3.9).Using Wick's decomposition of the operators, we obtain Here and henceforth we denote N − (q, t) = n(q, t), N + (q, t) = 1+n(q, t).The symbol "•" means the product of matrices, "Tr" denotes the trace over zone indices, h.c.stands for the hermitian conjugated term and I is the unit matrix.The correlation contributions to the collision integrals come from the last terms in equations (3.8) and (3.9).Some algebra (see Appendix B) gives for the corresponding contribution to the carrier collision integral where the matrix B(k, k − q, t) is defined as Explicit expressions for the Lagrange multipliers are derived in Appendix B. The correlation term in the phonon collision integral has a similar structure: Note that the kinetic equations for the distribution functions are incomplete in themselves.They must be supplemented by an equation of motion for the mean interaction energy Ĥint t or, equivalently, by an equation for the quasi-temperature T (t) = 1/β(t).To derive the equation for the quasi-temperature, we shall proceed in the similar manner as in [9].Differentiating the self-consistency conditions (3.3) with respect to time and using the von Neumann equation (2.11) for ρ(t), a little algebra leads to The last equation can be converted into the equation for β(t), taking into account that the relevant statistical operator (3.1) depends on time through the Lagrange multipliers and then eliminating ∂Λ αβ (k, t)/∂t, ∂Λ ph (q, t)/∂t.Details are described in [9] (see Appendix B) and we do not repeat the algebra here.In the case of a weak electron-phonon interaction, the resulting equation for the inverse quasi-temperature reads where C(t) = Ĥint , Ĥint t is the "energy-energy" correlation function in the relevant ensemble.We use the conventional notation for the correlation function: Here the averages are calculated with the "truncated" statistical operator , we obtain it in terms of the nonequilibrium distribution functions: Now, since the explicit expressions for the collision integrals are derived, the kinetic equations (2.12) and the evolution equation (3.20) for the quasi-temperature provide a closed description of the ultrafast dynamics in a laser excited semiconductor.In the next sections we shall discuss some important properties of these equations.

Energy balance equation
The energy balance equation follows directly from equations (3.19) and can be written in the form where ε kin (t), ε ph (t), and ε int (t) are, respectively, the mean kinetic energy of the carriers, the mean energy of phonons, and the mean interaction energy: The polarization current is given by If one makes any approximation for kinetic kernels (for instance, expansion in series in external field amplitude), balance equation (4.1) will be broken automatically.
In such a case, one can only speak about balance equation up to a certain order in E 0 .On the other hand, at the relaxation stage, the external laser field E(t) is zero and, consequently, equation (4.1) becomes the energy conservation law for the carrier-phonon system.It should be emphasized that the balance equation (4.1) is valid for arbitrary collision integrals.This is an important advantage of our approach over other schemes where the evolution equation for the interaction energy is not completely consistent with kinetic equations for carriers and phonons.Using the collision integrals (3.15)-(3.18)derived in non-Markovian Born approximation, one can obtain an "explicit" expression for the mean interaction energy in terms of the single-particle distribution functions and the quasi-temperature.To do this, we multiply the first equation (2.12) by the single-particle energy matrix M(k, t) and the second equation by ω LO .Then, using the above formulae for the collision integrals, the definition of the time-ordered exponents (3.13) and their symmetry properties (see Appendix A), we arrive at the balance equation (4.1) with the mean interaction energy where ε int (t 0 ) is the initial value of the mean interaction energy, and Previously the expression for ε int (t) analogous to (4.4) was derived in [9] for the special case of a one-zone semiconductor in the absence of the external field.

Relaxation stage of the evolution and the Markovian limit
In general, equations (2.12) and (3.20) can be solved only numerically because the collision integrals have a rather complicated matrix structure.In this section we consider the problem of relaxation towards equilibrium and make the following approximations: (i) The collision integrals for carriers and phonons are taken in the second non-Markovian Born approximation (see section 3); (ii) Off-diagonal elements are neglected in the kinetic kernels which corresponds to their field-free form; (iii) Polarization, being the off-diagonal distribution function, is supposed to relax more rapidly than the electron, hole and phonon distribution functions2 .
In our opinion, the last assumption is more reasonable than the phonon bath approximation [2] where the phonon collision integral is taken to be zero, while the carriers distribution functions and the polarization are essentially nonequilibrium.With the above approximations, the two-band problem reduces to the effective one-zone problems for the electrons and holes distribution functions coupled by the evolution of the quasi-temperature and the phonon distribution.After some algebraic manipulations with the collision integrals cited in section 3, the kinetic equations at the relaxation stage can be written as ∂n(q, t) ∂t = 2D 2 (q) where we have introduced the designations: . (5.4) Within the same approximations, the evolution equation (3.20) for the quasi-temperature takes the form with the following expression for C(t): Since {(x − 1)/ ln x} > 0, we see that C(t) > 0. It is instructive to show that the equilibrium carrier and phonon distribution functions are stationary solutions of equations (5.1) and (5.2).First we define the auxiliary quasiequilibrium distribution functions where the nonequilibrium chemical potential of the carriers, µ(t), is determined by the condition of electroneutrality It can easily be verified that where the K-function in the right-hand side is obtained from the (5.4) if the nonequilibrium distribution functions are replaced by the quasiequilibrium functions (5.7).Relation (5.9) may tell something useful about the properties of the evolution equations (5.1), (5.2), and (5.5).We see that the right-hand sides of these equations are so constructed that the collision and correlation contributions exactly cancel each other if f α (k, t ) = f 0 α (k, t ) and n(q, t ) = n 0 (t ).In particular, if the evolution of the system starts at t 0 with thermal equilibrium, the distribution functions and the quasi-temperature will remain unchanged for all times t > t 0 .It must be emphasized that the non-Markovian evolution equations without correlation terms do not have an equilibrium solution!Let us dwell briefly on the transition to the Markovian limit in the evolution equations.The Markovian approximation implies that the collision integrals in kinetic equations (5.1), (5.2), and the right-hand side of equation (5.5) depend on the distribution functions and the quasi-temperature taken at time t.The Markovian approximation is inadequate for the initial stage of evolution and, strictly speaking, describes only the limiting long-time regime where ∂f α (k, t)/∂t → 0, ∂n(q, t)/∂t → 0, and ∂β(t)/∂t → 0. On passage to the Markovian approximation in kinetic equations (5.1), (5.2), and in equation (5.5) for the quasi-temperature, we put f α (k, t ) ≈ f α (k, t), n(q, t ) ≈ n(q, t), and then take the limit t − t 0 → ∞.The remainder integrals over t are calculated in the standard way: (5.10) We see immediately that in the Markovian limit the correlational components of the collision integrals (the terms including β(t)), as well as the right-hand side of (5.5), become zero.Indeed, in all these expressions we obtain a construction like Ωδ(Ω) ≡ 0. Nevertheless, the "collision" components of I α , I ph remain finite and have the form of the Uehling-Uhlenbeck electron-phonon collision integrals (see [9,12]).Thus, the correlation terms do not appear in the Markovian collision integrals, but this is true only in the second Born approximation.In higher-order approximations, such terms do not vanish in the Markovian limit (see, e.g.[12]).

Inclusion of quasiparticle damping
An important point is that the collision integrals (3.15)-(3.18)lead to non-Markovian kinetic equations with infinite memory depth.This is a consequence of the fact that the evolution operator (3.7) appearing in the nonequilibrium statistical operator (3.6) corresponds to free particle approximation in the microscopic dynamics.It is clear, however, that kernels in kinetic equations must decay in time due to quasiparticle damping.That the quasiparticle damping plays an important part in non-Markovian kinetics of ultrafast processes is now well understood (see, e.g., [2,6] for a detailed discussion).The simplest improvement is to add a phenomenological small damping into the kernels (3.13), and try to determine γ as a fit parameter 3 .In more sophisticated approaches, the kinetic kernels can be calculated, in principle, from some fundamental quantities.
For example, in the Green's function formalism [2,6] the kinetic kernels are expressed in terms of retarded and advanced Green's functions or, equivalently, in terms of the spectral function.The choice (6.1) corresponds to a Lorentzian spectral function.Unfortunately, Lorentzian spectral functions do not conserve total energy, nor do they yield the correct asymptotic distribution functions at the relaxation stage.Although many attempts have been made to introduce damped kinetic kernel systematically [3,6,[13][14][15][16][17], the problem is still essentially unsolved, and the investigation of improved spectral functions is a field of ongoing research.We would like to emphasize that in our approach the specific form of the quasiparticle damping is not so important as in other available methods.As mentioned, if the mean interaction energy is treated as an independent state parameter, then the correct energy balance is valid for any collision integrals (see section 4).Different approximations for collision integrals I αβ and I ph result in different equations for the quasi-temperature [see (3.20)] and, consequently, in different approximations for the mean interaction energy.One may say that in our approach the energy balance is used as a constraint on the kinetic kernels and the mean interaction energy.It is thus hoped that even the simple choice of damped kinetic kernels, equation (6.1), with a fit parameter γ can yield the correct behavior of the distribution functions.

Numerical results
In our previous paper [18] we investigated the dynamics of a two-zone semiconductor during the laser pulse irradiation and the relaxation toward equilibrium with an additional assumption that the polarization relaxes much faster than the carrier and phonon distribution functions.Here we present some results for the relaxation stage of the evolution.It is assumed that the evolution of the system starts with an initially nonequilibrium state and the initial polarization is negligible 4 .As in the paper [18], we have chosen the following typical energies and time scales: the energy of longitudinal optical phonons ω LO =36 meV, the Rabi energy ω R =26.3 meV, the laser pulse duration τ pulse =50 fs, the width of the band gap E g =1.8 eV, the excess energy ∆ = ω − E g =60 meV.The dimensionless constant of the electron-phonon interaction was 0.069; the mass ratio of holes to electrons was 6.86 and the initial temperature of the phonon bath was 164 K.We use the values ω LO and τ pulse as the characteristic energy and the time scale to obtain dimensionless quantities.To solve the kinetic equations (5.1)-(5.2) and the equation for quasi-temperature (5.5), we use the factorization of the kinetic kernels.Details of the transition from integrodifferential equations to an extended system of ordinary differential equations are presented in Appendix C.This extended set is solved by applying the fifth order Runge-Kutta integration with automatic step-size control.We chose a dense enough energy grid with N = 60 points that corresponds to ∼ 10 6 ordinary differential equations.Figures 1 and 2 show the energy dependence of the electron distribution functions at fixed times in the case of γ = 0.It is seen that the neglect of the nonequilibrium correlations (figure 1) leads to nonphysical behavior of the electron distribution function at long times: f c (ε, t) exceeds unity and is far away from its quasiequilibrium asymptote.Correlation term I (corr) improves the situation greatly (figure 2) though does not ensure complete relaxation to the equilibrium (see location of the lines with stars and half-filled pentagons at low and high energies).In other words, it means that the second Born approximation is not sufficient to establish equilibrium at long times.Therefore we have to introduce a nonzero γ, which formally means going beyond the second Born approximation.In figure 3 we presented long time behavior of the electron distribution function with taking into consideration both nonequilibrium correlations and nonzero quasiparticle damping 5 .One can see complete relaxation of the distribution function to its asymptotic value (a line with triangles at t * = 40τ pulse ).Thus, our numerical results are in accord with relaxation dynamics of [2,3]; moreover, our approach enables one to analyze the asymptotics of phonon distribution functions (see figure 5) by means of quasi-temperature.Figure 4 shows the time dependence of the total phonon number at the relaxation stage.On the left panel in figure 4 one can that collision integrals in the 2-nd Born approximation without any relaxation mechanisms (see explanation does not ensure relaxation of the total phonon number towards equilibrium.However, once we have introduced an additional term I add ph (q, t) = (n 0 (t) − n(q, t))/T ph (it corresponds to the phenomenological phonon lifetime) to the phonon collision integral in a form of a simple τ -approximation (right panel), we obtained a relaxation of N ph (t) = q n(q, t) to its asymptotic value.Moreover, as seen in figure 4b, the correlation term in the phonon collision integral results in a smooth relaxation (solid line), while the neglect of the running nonequilibrium correlations leads to non-motivated flexure of a dashed line at long times.However, a question about an appropriate choice of T ph arises.In our calculations, the phonon lifetime was treated as a fit parameter to ensure relaxation of n(q, t) to equilibrium.
Figure 5 illustrates the long time behavior of the phonon distribution functions.It is seen that n(q, t) is close to its quasiequilibrium value n ph (ε, t)| t τ pulse = n 0 ph (t) = 1/[exp( ω LO /T (t)) − 1] shown by the dashed line.However, there is a deviation of the stationary distribution of the phonons from its quasiequilibrium value.Whereas this deviation is only about a few percent for zero γ and a relatively large phonon lifetime T ph = 7τ pulse , it reaches much more pronounced values when a quasiparticle damping increases to 1/10 and the phonon relaxation time decreases to 3τ pulse .In our opinion, this is connected with a lack of information about the dynamics of quasiparticle damping and could serve as a manifestation of the necessity of considering nonequilibrium correlations and the formation of quasiparticle behavior on equal terms.

Conclusions
In this paper we have derived a closed set of evolution equations for a two-zone semiconductor excited by an ultrashort laser pulse, taking into account nonequilibrium correlations in the system.These evolution equations posses two basic properties: they have an equilibrium solution and energy balance.It was shown that the introduction of nonequilibrium correlations is of great importance.The corresponding terms in the collision integrals ensure the existence of the Markovian limit and equilibration of the system.However, the numerical of the relaxation stage show that the 2-nd Born approximation is not sufficient for a proper description of the at time scales.This is not strange, because the effect of higher order interactions is known to be essential at the stage of the formation of quasiparticle behavior [3,19,20].We found that even a rather small constant quasiparticle damping can stabilize the solutions of the kinetic equations and ensure relaxation of the distribution functions to their long time asymptotics.In the general case it is necessary to compose the time evolution equation for guasiparticle damping, which supplements the system of kinetic equations together with the dynamic equation for quasi-temperature.To describe the dynamics of the system during the laser excitation, as well as at the relaxation stage, one has to consider a non-simplified version of the total collision integrals for carriers, phonons and polarization including the field contributions and quasiparticle damping.We plan to perform such an extended analysis in future publications.
strong time-dependent field, these equations cannot be solved analytically.To find explicit expressions for the kinetic kernels, we neglect the field-dependent terms in the equations for diagonal elements.Then we find that Substituting these expressions into equations (3.14) for the off-diagonal elements and assuming that the field amplitude varies with time much than the cosine factor in (2.8), which allows us to perform integration over t explicitly, we obtain the following form of the field-dependent kinetic kernels: where δ − k = e ck − e vk − ω denotes detuning of the system6 .Let us make some conclusions from the expressions (A.1)-(A.3).On the one hand, the kinetic kernels possess the symmetric properties similar to those of retarded and advanced Green's functions.On the other hand, the series in a field amplitude turn out to be factorizable in t, t : The Lagrange multipliers Λ αβ (k, t) and Λ ph (q, t) can be eliminated with the aid of the self-consistency conditions (3.3).Evaluation of the averages with the truncated relevant statistical operator (3.22) gives [exp{Λ(k, t)} + I ] −1 αβ = f αβ (k, t), exp{Λ ph (q, t)} − 1 −1 = n(q, t).(B.5) It is immediately obvious that Λ ph (q, t) = ln 1 + n(q, t) n(q, t) .(B.6) with the initial conditions G i (e k , e k , e k , t 0 ) = 0. Thus, we arrive at the system of ordinary differential equations where the energies e k play the role of parameters.Then, passing from the integrals over e k , e k to the sums with the corresponding Gaussian weights, we obtain a system of dim = 5N + 12N 3 + 1 equations.Here 5N equations are for the one-particle distribution functions (for carriers, phonons and polarization), 12N 3 equations describe the dynamics of the auxiliary kernels (C.2), and the last equation determines the quasi-temperature T (t).
Taking the energy grid with N = 60 points, we obtain the system of about 10 6 equations which are to be solved numerically with given initial conditions.

Figure 1 .
Figure 1.The electron distribution function f c (ε, t) at fixed times.Nonequilibrium correlations are not taken into account.

Figure 2 .
Figure 2. The electron distribution function f c (ε, t) at fixed times.Nonequilibrium correlations are taken into account.

Figure 3 .
Figure 3.The electron distribution function f c (ε, t) at fixed times.Nonequilibrium correlations are taken into account and nonzero quasiparticle damping γ = 0.1 is introduced.