text:
# On dispersive energy transport and relaxation in the hopping regime

## I Introduction

In recent years much attention has been devoted to the study of relaxation processes of non-equillibrium charge carriers in strongly localized systems, where transport proceeds via phonon-assisted hopping, like photoexited charge carriers in band tails (see e.g. -) and Anderson insulators (see e.g.-). In such systems particular small relaxation speeds are observed. Often the smallness of the relaxation speed is attributed to interaction effects. However, also in strongly localized non-interacting electron systems long lasting relaxation processes are known to be the rule, and not the exception .

The theoretical investigation of such relaxation processes is notoriously difficult, since the system is strongly disordered and always in a transient state. In most problems one is interested in time scales which are large as compared to the time a single hop needs. For such time scales all quantities depend strongly on frequency, also for very low frequencies, so that the consideration of dispersive transport is vital. On the other hand in most problems of interest both spatial and energetical disorder exists. Due to the latter fact the transitions are inelastic. The inelastic character of the transitions, the relaxation, leads to a flow of energy from the electron system to the phonon system. Due to disorder this transport of energy is also dispersive. Therefore, dispersive energy transport and relaxation are intimately connected. The investigation of the relaxation process requires both the consideration of dispersive particle transport and the consideration of dispersive energy transport.

The intricate physical situation manifests also in the equations, which have to be considered. Since the transport is inelastic already in the simplest approximation the investigation requires to find solutions to integral equations (see e.g.\- ). Due to the fact that the particle moves in an energy dependent density of states the kernel of these integral equations does not depend on the difference between the site energies only. The traditional method to handle the situation, the percolation theory, is not applicable here. The effective-medium methods by Movaghar and coworkers give, as pointed out by Movaghar and coworkers, wrong results for systems at low temperature  . Furthermore, beyond the Markovian approximation the derivation of these integral equations represents itself an intricate problem. To our knowledge, beside the attempts by Movaghar and coworkers, numerical investigations (see e.g., ), and physical intuitive considerations (see e.g. ,), mainly Markovian equations have been used (see e.g. ).

It is the aim of the present paper to fill this gap and to provide a formalism that can be used for studying relaxation phenoma of strongly localized, charge carriers far from equillibrium in the hopping regime taking into account both dispersive particle transport and dispersive energy transport. To this end we focus, for the sake of definitness, on relaxation of photoxited, non-interacting, non-equillibrium charge carrires in band tails of, e.g. amorphous semiconductors, like amorphous Si:H. On the first glance, this problem seems to be rather special. If the number of charge carriers exited is small Fermi-correlation is negligible, so that one has to cope with linear rate equations. The linearity of the transport equations is, of course, the basic ingredient in solving the problem. A closer look on the problems of interest reveals, however, that most of the problems can be formulated in this way. Clearly, how to achieve linearization depends on the experiment chosen, but, provided the number of charge carrires exited is small, always the smallness of the particle number can be invoked. In this case, the number of particles at a site is certainly not small as compared to its equillibrium value, but small as compared to unity. After linearization the structure is, in principle, quite similar as for the case considered here.

Below we present the derivation of a framework for the consideration of relaxation phenomena due to phonon-assisted hopping at zero temperature. We focus on the limit of strong localization, where dispersive effects are expected to be most pronounced. Here the disorder manifests in a strong dependence of the transport coefficients on frequency for very low frequencies. Consequently, in this regime the diffusion propagator can not be calculated from Markovian transport equations, like those used e.g. in Ref.. To simplify our integral equations we use the concept of quasi-elasticity, the particle changes its energy only by a small amount by one hop.

The method is applied to the study of dispersive energy transport. We find that always two tendencies are present, if only the density of states not decreases with increasing energy. First, there is a widening of the packet due to the statistical spreading, the dispersive energy transport. Second, there is narrowing of the packet due to decrease of the density of states with decreasing energy. The overall evolution is determined by the interplay of these two tendencies. We have studied both tendencies in two limiting cases, for a particle moving in a constant density of states and for a particle moving in an exponential density of states. In the first case the impact of the density of states decrease is absent, so that only the statistical spreading is present. Here the packet evolves into a Gaussian packet moving with constant drift velocity in energy space. The packet width increases with time as $`\sqrt{t}`$. The other result is obtained for exponential density of states. Here both tendencies are present. We find, that the velocity of the packet is strongly slowing down with time. In that case the mean square deviation of energy eventually becomes time-independent. Consequently, the motion of the packet in energy space is “soliton-like”. The concrete results on the diffusion propagator, its time dependence and its moments for exponential densities of states are of relevance for photoluminiscense experiments on amorphous Si:H.

## II Basic equations

We consider photoexited, localized charge carriers in band tails at zero temperture. After exitation the charge carriers lower their energy by phonon-assisted hops between localized states. Since $`T=0`$, only hops from higher to lower energy occur. In this situation the charge carriers are in strong non-equilibrium. We assume, that the number of excited charge carriers is small, and their energies are far from the Fermi level, so that it is very unlikely that an electron jumps to a site already occupied. Consequently, we can neglect Fermi-correlation. In this case the electron transport can be described by the rate equation .

$$\frac{dp_m}{dt}=\underset{m^{}}{}\left(p_m^{}W_{m^{}m}p_mW_{mm^{}}\right).$$
(1)
In calculating the transition probabilities we assume, that the electron-phonon coupling strength is weak, so that only one-phonon processes have to be taken into account. Then, at zero temperature, the transition probabilities are given by

$$W_{m^{}m}=\mathrm{\Theta }\left(\omega \epsilon _m^{}+\epsilon _m\right)\mathrm{\Theta }\left(\epsilon _m^{}\epsilon _m\right)W\left(\left|𝐑_{mm^{}}\right|\right),$$
(2)
where

$$W\left(\left|𝐑_{mm^{}}\right|\right)=\nu e^{2\alpha \left|𝐑_{mm^{}}\right|}.$$
(3)
Here $`\alpha `$ is the inverse of localized state radius, and $`\nu `$ is the phonon frequency. The energy $`\omega `$ is the upper bound for the energy transfer by one hop. Note, that in the materials of interest $`\omega `$ can be much smaller as the Debye-frequency, since not all phonons can interact with localized electrons equally well. Short wave-length phonons are ineffective since the electron-phonon coupling constant tends to zero for momenta $`q`$ with $`q/2\alpha 1`$. Therefore, the effective upper phonon momentum is of the order $`2\alpha `$, and not of the order of the inverse lattice constant of the host material. Furthermore, in disordered systems the high-energetic phonons are localized, and need not contribute to transport by necessity.

The first step function in front of the transition probabilities restricts the transitions to transitions between sites separated by at most $`\omega `$ in energy space. Thus, it decreases the energy relaxation speed. In impurity conduction, and in nearly all papers on relaxation of charge carriers in band tails, this step function is usually replaced by unity. In impurity conduction this is quite reasonable, since it is assumed that hops are restricted within narrow stripe near the Fermi level, which is small compared to Debye energy. In the band-tail problem, however, we can see no reason for neglecting it in advance.

To calculate the transport quantities of interest, we have to calculate the diffusion propagator. In order to render the analytical calculations feasible we introduce continuous coordinates. The change of representation is defined by:

$$n(\rho )=\underset{m}{}p_m\delta \left(\rho \rho _m\right),$$
(4)
where $`\rho =(𝐑,\epsilon )`$ and $`\rho _m=(𝐑_m,\epsilon _m)`$. In this representation the rate equation (1) takes the form:

$$\frac{dn\left(\rho \right)}{dt}=𝑑\rho ^{}n\left(\rho ^{}\right)V(\rho ^{},\rho ),$$
(5)
where $`V`$ is determined by the equations:

$$V(\rho ^{},\rho )=𝑑\rho _1\eta \left(\rho _1\right)w_{\rho _1}(\rho ^{},\rho ),$$
(6)
$$w_{\rho _1}(\rho ^{},\rho )=W(\rho ^{},\rho _1)\left[\delta \left(\rho _1\rho \right)\delta \left(\rho ^{}\rho \right)\right],$$
(7)
$$\eta \left(\rho \right)=\underset{m}{}\delta \left(\rho \rho _m\right).$$
(8)
The Laplace-transformed equation is given by:

$$sn\left(\rho \right)n_o\left(\rho \right)=𝑑\rho ^{}n\left(\rho ^{}\right)V(\rho ^{},\rho ),$$
(9)
where $`n_o(\rho )=n(\rho ,t=0)`$ is the initial condition. We will assume that $`p_m(t=0)`$ is a function $`p_0(𝐑_m,\epsilon _m)`$, so that $`n_o(\rho )=p_0(\rho )\eta (\rho )`$.

Equation (9) can be solved using the Green’s function method. The solution is given by:

$$n\left(\rho \right)=𝑑\rho ^{}n_o\left(\rho ^{}\right)\mathrm{\Phi }(\rho ^{},\rho ).$$
(10)
The Green’s function satisfies the equation:

$$s\mathrm{\Phi }(\rho ^{},\rho )𝑑\rho _1V(\rho ^{},\rho _1)\mathrm{\Phi }(\rho _1,\rho )=\delta \left(\rho ^{}\rho \right).$$
(11)
Note, that, due to probability conservation, the Green’s function $`\mathrm{\Phi }`$ and the probability $`w_{\stackrel{~}{\rho }}`$ satisfy the relations:

$$𝑑\rho \mathrm{\Phi }(\rho ^{},\rho )=\frac{1}{s},$$
(12)
$$𝑑\rho w_{\stackrel{~}{\rho }}(\rho ^{},\rho )=0.$$
(13)
## III Configuration Average

In order to calculate the configuration average, we assume, that the sites are distributed homogeneously in space. The distribution of site energies is supposed to be given by some distribution function $`p(\{\epsilon _i\})`$. Accordingly, the average of any quantity, depending on the energy and positions of sites, is given by:

$$A=\mathrm{\Pi }_m\frac{d𝐑_m}{𝒱}𝑑\epsilon _mp\left(\left\{\epsilon _m\right\}\right)A\left(\{𝐑_m,\epsilon _m\}\right),$$
(14)
where $`𝒱`$ is the volume of the system. The application of the averaging procedure to the structural factor $`\eta `$ serves, in particular, as a definition for the density of states, i.e.

$$𝒩\left(\epsilon \right)=\eta \left(\rho \right).$$
(15)
Products of the structural factor $`\eta `$ are averaged according to:

$$\eta \left(\rho _1\right)\mathrm{}\eta \left(\rho _n\right)=𝒩\left(\rho _1\right)\delta \left(\rho _1\rho _2\right)\mathrm{}\delta \left(\rho _{n1}\rho _n\right).$$
(16)
Using these definitions, the configuration average of Eq.(10) can be calculated diagrammatically . The diagrammatic method leads to the following set of equations for the calculation of the configuration average $`n(\rho )`$ of the electron density $`n(\rho )`$ , (see Figs 14 for illustration):

$$n\left(\rho \right)=𝑑\rho _1𝑑\rho _2p_0\left(\rho _1\right)S(\rho _1,\rho _2)F(\rho _2,\rho ),$$
(17)
$$sF(\rho ^{},\rho )=\delta \left(\rho ^{}\rho \right)+𝑑\rho _1\mathrm{\Pi }(\rho ^{},\rho _1)F(\rho _1,\rho ),$$
(18)
$$\mathrm{\Pi }(\rho ^{},\rho )=𝑑\rho _1𝒩\left(\rho _1\right)\mathrm{\Pi }_{\rho _1}(\rho ^{},\rho ),$$
(19)
$`\mathrm{\Pi }_{\rho _1}(\rho ^{},\rho )=w_{\rho _1}(\rho ^{},\rho )`$ (20)
$`+`$ $`{\displaystyle 𝑑\rho _2𝑑\rho _3w_{\rho _1}(\rho ^{},\rho _2)F(\rho _2,\rho _3)\mathrm{\Pi }_{\rho _1}(\rho _3,\rho )},`$ (21)
$$S(\rho ^{},\rho )=𝒩\left(\epsilon ^{}\right)\left[\delta \left(\rho ^{}\rho \right)+𝑑\rho _1F(\rho ^{},\rho _1)\mathrm{\Pi }_\rho ^{}(\rho _1,\rho )\right].$$
(22)
Here $`F(\rho ^{},\rho )=\mathrm{\Phi }(\rho ^{},\rho )`$. A detailed derivation of the set of integral equations is given in Appendix A.

## IV Effective Medium

Given the system of integral equations (18)-(22), the main problem is to find an approximate self-consistent solution to it. The situation is quite similar to that of the calculation of the equilibrium conductivity in a disordered system . There an approximate solution of this system could be found by introducing a proper decomposition of the function $`F`$, the diffusion propagator, into short- and long-wavelength limit, according to:

$$F(\rho ^{},\rho )=f(s)C(\epsilon )\delta \left(\rho ^{}\rho \right)+\stackrel{~}{F}(\rho ^{},\rho ),$$
(23)
where $`f(s)`$ was a frequency dependent parameter that could be related to the critical hopping length $`R_c`$ via the equation:

$$f\nu =\mathrm{exp}\left(2\alpha R_c\right),$$
(24)
and $`C(\epsilon )`$ was an energy dependent function determined by the principle of detailed balance. This decompostion originated from the notion, that the integrals in the integral equations are governed by products of transition probabilities and diffusion propagators, and, for strongly localized electrons, the latter quantities are short ranged functions as compared to the transition probabilities. Using the decomposition (23), the effective medium approximation reduces to the repacement of $`F`$ by $`fC(\epsilon )\delta (\rho ^{}\rho )`$ in the calculation of the effective transition probability $`\mathrm{\Pi }`$ (Eq.(21)) and of the irreducible block $`S`$ (Eq.(22)).

Here in the band-tail problem we use the same philosophy. We first decompose the diffusion propagator in two parts, according to:

$$F(\rho ^{},\rho )=f(\epsilon ,s)\delta \left(\rho ^{}\rho \right)+\stackrel{~}{F}(\rho ^{},\rho ).$$
(25)
Then, to investigate the consequences of this decomposition, we insert Eq.(25) into Eq.(21). Performing a partial summation we obtain:

$`\mathrm{\Pi }_{\stackrel{~}{\rho }}(\rho ^{},\rho )=\stackrel{~}{w}_{\stackrel{~}{\rho }}(\rho ^{},\rho )`$ (26)
$`+`$ $`{\displaystyle 𝑑\rho _1𝑑\rho _2\stackrel{~}{w}_{\stackrel{~}{\rho }}(\rho ^{},\rho _1)\stackrel{~}{F}(\rho _1,\rho _2)\mathrm{\Pi }_{\stackrel{~}{\rho }}(\rho _2,\rho )},`$ (27)
where the renormalized transition probabilities $`\stackrel{~}{w}_{\stackrel{~}{\rho }}(\rho ^{},\rho )`$ are given by Eq.(7), with $`W`$ replaced by:

$$\stackrel{~}{W}(\rho ^{},\rho ;s)=\frac{W(\rho ^{},\rho )}{1+f(\epsilon ^{},s)W(\rho ^{},\rho )}.$$
(28)
Note that the renormalized transition probability depends now on $`s`$ via the function $`f(\epsilon ,s)`$.

At this point the introduction of the function $`f(\epsilon ,s)`$, the effective medium, still seems to be rather arbitrary. However, if we now choose the effective medium in such a way, that the integrals over $`\stackrel{~}{F}`$ vanish, the renormalized transition probabilities yield the exact solution of the diffusion problem. In that case the function $`S`$ turns into:

$$S(\rho ^{},\rho )=𝒩\left(\epsilon ^{}\right)\delta \left(\rho ^{}\rho \right),$$
(29)
so that Eq.(17) can be cast into the form:

$$n(\rho ;s)=𝑑\rho ^{}n_o\left(\rho ^{}\right)F(\rho ^{},\rho ).$$
(30)
Therefore the function $`F(\rho ^{},\rho )`$ can be identified with the diffusion propagator.

Note, that although the general philosophy is quite the same as in equilibrium, the situation is much more intricate here, — the principle of detailed balance is absent, since we are dealing with the situation at zero temperature, and so it can’t determine the energy dependence of the effective medium. Moreover, we expect the critical hopping length to depend somehow on the position of the electron in the tail. Thus, in contrast to equilibrium, $`f`$ is not only a frequency dependent parameter but also an energy dependent function, which has to be determined self-consistently.

## V Diffusion propagator in effective medium approximation

For the moment we put aside the question of the determination of the effective medium to elaborate further on the consequences of the renormalization of the transition probabilities. To this end we focus on the diffusion propagator.

The equation for the diffusion propagator is given by Eq.(14). In effective medium approximation, when calculating the irreducible part $`\mathrm{\Pi }`$, only the lowest order contribution to $`\mathrm{\Pi }`$ with respect to $`\stackrel{~}{F}`$ is taken into account, so that $`\mathrm{\Pi }_{\stackrel{~}{\rho }}=\stackrel{~}{w}_{\stackrel{~}{\rho }}`$. In that approximation the equation for the diffusion propagator in momentum representation reads:

$$sF\left(𝐪|\epsilon ^{},\epsilon \right)=\delta \left(\epsilon ^{}\epsilon \right)+𝑑\epsilon _1\left\{F\left(𝐪|\epsilon ^{},\epsilon _1\right)\stackrel{~}{W}\left(𝐪|\epsilon _1,\epsilon ;s\right)𝒩\left(\epsilon \right)F\left(𝐪|\epsilon ^{},\epsilon \right)\stackrel{~}{W}\left(\mathrm{𝟎}|\epsilon ,\epsilon _1;s\right)𝒩(\epsilon _1)\right\}.$$
(31)
Of course, as we don’t know how the effective medium looks like so far, and moreover, as the equation is a complicated integral equation, we can’t find a solution. The fact, that the effective medium is a function of energy, makes the problem much more complicated. Progress can only by achieved if we can find arguments to simplify the equation considerably. To this end we focus on the renormalized transition probability. According to the Eqs. (2) and (28), it is given by:

$$\stackrel{~}{W}\left(R|\epsilon ^{},\epsilon ;s\right)=\mathrm{\Theta }\left(\epsilon ^{}\epsilon \right)\mathrm{\Theta }\left(\omega \epsilon ^{}+\epsilon \right)\frac{W\left(R\right)}{1+f(\epsilon ^{},s)W(R)}=\mathrm{\Theta }\left(\epsilon ^{}\epsilon \right)\mathrm{\Theta }\left(\omega \epsilon ^{}+\epsilon \right)\stackrel{~}{W}\left(R|\epsilon ^{};s\right).$$
(32)
Owing to the step functions in front of the transition probability, the energy integrations are restricted to intervals of length $`\omega `$. Taking this fact into account the integral equation for the calculation of the diffusion propagator takes the form:

$$sF(𝐪|\epsilon ,\epsilon )=\delta (\epsilon ^{}\epsilon )+\underset{0}{\overset{\omega }{}}d\epsilon _1\{F(𝐪|\epsilon ^{},\epsilon +\epsilon _1)\stackrel{~}{W}(𝐪|\epsilon _1+\epsilon ;s)𝒩\left(\epsilon \right)F(𝐪|\epsilon ^{},\epsilon )\stackrel{~}{W}(0|\epsilon ;s)𝒩(\epsilon _1+\epsilon )\}.$$
(33)
Further we assume that the diffusion propagator, the transition probability $`\stackrel{~}{W}(R,\epsilon ;s)`$, and the density of states $`𝒩\left(\epsilon \right)`$ are slowly varying functions on intervals of length $`\omega `$, so that the integrand can be expanded with respect to $`\epsilon _1`$. Doing so we obtain:

$$sF(𝐪|\epsilon ,\epsilon )=\delta (\epsilon ^{}\epsilon )+\omega 𝒩\left(\epsilon \right)[\stackrel{~}{W}(𝐪|\epsilon ;s)\stackrel{~}{W}(\mathrm{𝟎}|\epsilon ;s)]F(𝐪|\epsilon ^{},\epsilon )+\frac{1}{2}\omega ^2\frac{}{\epsilon }\left[F(𝐪|\epsilon ^{},\epsilon )\stackrel{~}{W}(𝐪|\epsilon ;s)𝒩\left(\epsilon \right)\right].$$
(34)
We terminate the expansion after the first derivative with respect to energy. This term describes the biased motion of the particle from sites of higher to sites of lower energy. For finite temperatures one would have to replace $`F\left(𝐪|\epsilon ^{},\epsilon \right)`$ in the last term with $`F\left(𝐪|\epsilon ^{},\epsilon \right)+kTF\left(𝐪|\epsilon ^{},\epsilon \right)/\epsilon `$ which would describe the energy diffusion current. At zero temperature there is no input of energy from the phonon-system into the particle. Therefore, we expect the additional term to be rather small, and can be neglected.

To investigate the transport properties, we restrict our attention with the long-wavelength limit. In that case the elastic contribution yields the diffusion coefficient for particle diffusion. The term, containing the derivative with respect to energy is finite for $`q0`$. If we are interested only in the long-wavelength limit, we can put here $`q=0`$, since the remaining terms are of higher order with respect to ’$`q`$’ and ’$`\omega /\epsilon `$’. Since a nonzero momentum in this term couples particle diffusion to energy transport, this approximation corresponds to a decoupling of these two processes. Then, in the long-wavelength limit, we obtain:

$`sF\left(𝐪|\epsilon ^{},\epsilon \right)`$ $`=`$ $`\delta \left(\epsilon ^{}\epsilon \right)D(\epsilon ,s)q^2F\left(𝐪|\epsilon ^{},\epsilon \right)`$ (35)
$`+`$ $`{\displaystyle \frac{}{\epsilon }}\left[F\left(𝐪|\epsilon ^{},\epsilon \right)v(\epsilon ,s)\right],`$ (36)
where:

$$D(\epsilon ,s)=\frac{1}{2}\omega _q^2\stackrel{~}{W}\left(q|\epsilon ;s\right)|_{𝐪=\mathrm{𝟎}},$$
(37)
and:

$$v(\epsilon ,s)=\frac{1}{2}\omega ^2𝒩\left(\epsilon \right)\stackrel{~}{W}\left(0|\epsilon ;s\right).$$
(38)
Eq.(36) is easily solved. It’s solution is:

$$F\left(𝐪|\epsilon ^{},\epsilon \right)=\frac{\mathrm{\Theta }\left(\epsilon ^{}\epsilon \right)}{v(\epsilon ;s)}\mathrm{exp}\left[\underset{\epsilon }{\overset{\epsilon ^{}}{}}𝑑\epsilon _1\frac{s+D(\epsilon _1;s)q^2}{v(\epsilon _1;s)}\right].$$
(39)
Explicit expressions for the transport coefficients can be found in Appendix B.

Now, at this stage, the validity of our quasi-elastic approximation requires, that the second derivative terms are small as compared to terms with first derivatives with respect to energy. This requirement imposes the following restrictions on the transport coefficients $`v`$ and $`D`$, the frequency $`s`$ and the momentum $`q`$:

$$\omega \left|\frac{1}{v}\frac{dv}{d\epsilon }\right|1,$$
(40)
$$\omega \left|\frac{s}{v}\right|1,$$
(41)
$$q^2\omega \frac{D}{v}1.$$
(42)
The applicability of the quasi-elastic approximation was also discussed in Ref.. There it was concluded that this approximation should be inapplicable. To substantiate this statement numerical calculations were invoked. However, in interpreting these data it has to be taken into account, that they have been obtained using a model, neglecting a weighting of the transition probabilities according to the number of phonons emitted, so that hops between nearly isoenergetic sites were treated as likely as hops from the very top to the very bottom of the tail. Consequently the discussion in Ref. applies only to system with sufficient strong electron-phonon interaction, but not to systems with weak electron-phonon interaction. For extraordinary deep hops to contribute to the diffusion propagator they should be characteristic for an ensemble of electrons. This is, however, not expected and in the experimental data, e.g. on amorphous Si:H , not observed. For these reasons the discussion in Ref. does not apply to our model.

## VI Physical quantities and interpretation

Before establishing a self-consistency equation, we first elaborate further on the physical content of our diffusion equation. To this end let us first have a closer look on the coefficients $`D(\epsilon ,s)`$ and $`v(\epsilon ,s)`$. Imagine, that we have a particle initially located at ($`𝐑_0,\epsilon _0`$). Then the initial condition is $`n_o(𝐑,\epsilon )=\delta (𝐑𝐑_0)\delta (\epsilon \epsilon _0)`$. According to our approximation, the motion of the charge carrier is composed of two contributions: particle diffusion between isoenergetic sites, and relaxation in energy space. Characteristics of these two processes are the mean square displacement and the energy relaxation speed. It turns out, that both can be calculated from the function:

$`P_L(\epsilon _0,\epsilon ;s)={\displaystyle \frac{\mathrm{\Theta }\left(\epsilon _0\epsilon \right)}{s}}{\displaystyle \frac{}{\epsilon }}\mathrm{exp}\left[s{\displaystyle \underset{\epsilon }{\overset{\epsilon _0}{}}}{\displaystyle \frac{d\epsilon _1}{v(\epsilon _1,s)}}\right]`$ (43)
$`={\displaystyle \frac{\mathrm{\Theta }\left(\epsilon _0\epsilon \right)}{v(\epsilon ,s)}}\mathrm{exp}\left[s{\displaystyle \underset{\epsilon }{\overset{\epsilon _0}{}}}{\displaystyle \frac{d\epsilon _1}{v(\epsilon _1,s)}}\right],`$ (44)
which is just $`F(q=0|\epsilon _0,\epsilon )`$.

We define the mean energy of the particle by:

$$E\left(s\right)=𝑑\rho ^{}\epsilon ^{}n(\rho ^{},s).$$
(45)
Using our diffusion propagator, this equation can be rewritten in the form:

$$E\left(s\right)=𝑑\epsilon ^{}P_L(\epsilon _0,\epsilon ^{};s)\epsilon ^{},$$
(46)
and from the diffusion equation we obtain after integration by parts:

$$sE\left(s\right)\epsilon _0=𝑑\epsilon P_L(\epsilon _0,\epsilon ;s)v(\epsilon ,s).$$
(47)
Therefore, in general, the time dependence of the mean energy is given by complicated integrals. These integrals simplify considerably in two limiting cases: in the absence of dispersive energy transport, and for energy independent $`v(\epsilon ,s)`$. In the absence of dispersive energy transport, i.e. in the Markovian limit in which the transport coefficients are independent of $`s`$, the integrals can readly be calculated in time representation. In the latter situation we simply obtain:

$$P_L(\epsilon _0,\epsilon ;t)=\mathrm{\Theta }\left(\epsilon _0\epsilon \right)\delta \left(\epsilon _m\left(t\right)\epsilon \right),$$
(48)
where $`\epsilon _m(t)`$ is defined by:

$$t=\underset{\epsilon _m}{\overset{\epsilon _0}{}}\frac{d\epsilon _1}{v\left(\epsilon _1\right)}.$$
(49)
Therefore we obtain:

$$\frac{dE\left(t\right)}{dt}=\frac{d\epsilon _m\left(t\right)}{dt}=v\left(\epsilon _m\left(t\right)\right).$$
(50)
Consequently, $`v`$ is the velocity of energy relaxation, and $`\epsilon _m(t)`$ is the instantenous position of the particle in energy space. In general, however, in disordered systems the coefficient $`v(\epsilon ,s)`$ depends on $`s`$, so that energy transport is dispersive. In that case the integral in Eq.(47) can only be calculated easily if $`v`$ is independent of energy. Then, owing to probability conservation,

$$\left(\frac{dE}{dt}\right)\left(s\right)=\frac{v\left(s\right)}{s}$$
(51)
is obtained.

The mean square displacement, defined by:

$$R^2\left(t\right)=d\rho ^{}R^{}{}_{}{}^{2}n\left(\rho ^{}\right),$$
(52)
can be obtained from similar arguments. In Fourier representation this equation can be written in the form:

$$\left(\frac{dR^2}{dt}\right)\left(s\right)=s\mathrm{\Delta }_q|_{q=0}𝑑\epsilon ^{}F\left(q|\epsilon _0,\epsilon ^{}\right).$$
(53)
The derivative of the diffusion propagator can be calculated using the diffusion equation. Doing so, we obtain:

$$\left(\frac{dR^2}{dt}\right)\left(s\right)=2d𝑑\epsilon P_L(\epsilon _0,\epsilon ;s)D(\epsilon ,s).$$
(54)
Again Eq.(54) simplifies only in two special cases, in the Markovian limit and for energy independent diffusion coefficients. Whereas in the first case

$$\frac{dR^2}{dt}\left(t\right)=2dD\left(\epsilon _m\left(t\right)\right)$$
(55)
is obtained, we have

$$\left(\frac{dR^2}{dt}\right)\left(s\right)=2d\frac{D\left(s\right)}{s}$$
(56)
in the latter situation. Below we shall see, that energy independent transport coefficients are obtained for constant density of states only.

## VII The self-consistency equation

So far we have only investigated the consequences of the renormalization. In order to complete the approximation scheme, we still have to calculate the effective medium itself. Of course, we can not calculate the effective medium exactly, since this would amount to find an exact solution to the diffusion problem. Rather we shall try to calculate $`f(\epsilon ,s)`$ self-consistently.

The transport coefficients are properties of the Green’s function $`F`$, the diffusion propagator. Thus, in order to find an equation for $`f(\epsilon ,s)`$ we should relate the transport coefficients to $`\mathrm{\Pi }`$, the irreducible part of the diffusion propagator. The diffusion coefficient comprises only elastic contributions, however, in the relaxation problem a proper description of inelastic processes is vital. Moreover, as shown in the previuos section, the characteristics for particle and energy transport can already be calculated from the function $`P_L`$, the $`q0`$ limit of the diffusion propagator. Therefore, in establishing a self-consistency equation we should focus on those characteristics that are important in that limit, that is the energy relaxation speed $`v`$.

The equation for the diffusion propagator for $`q=0`$ is given by:

$`sP_L(\epsilon ^{},\epsilon ;s)=\delta \left(\epsilon ^{}\epsilon \right)+P_L(\epsilon ^{},\epsilon ,s){\displaystyle 𝑑\epsilon _1\mathrm{\Pi }(\epsilon _1,\epsilon ;s)}`$ (57)
$`+{\displaystyle \frac{P_L(\epsilon ^{},\epsilon ;s)}{\epsilon }}{\displaystyle 𝑑\epsilon _1\left(\epsilon _1\epsilon \right)\mathrm{\Pi }(\epsilon _1,\epsilon ;s)}.`$ (58)
If we compare Eq.(58) with Eq.(36) we deduce that

$$v(\epsilon ,s)=𝑑\epsilon _1\left(\epsilon _1\epsilon \right)\mathrm{\Pi }(\epsilon _1,\epsilon ;s).$$
(59)
Now we decompose $`v`$, as defined in Eq.(59), into two parts, one part that contains only the effective medium approximation of $`\mathrm{\Pi }`$:

$$v(\epsilon ,s)=𝑑\epsilon _1\left(\epsilon _1\epsilon \right)\mathrm{\Pi }|_{EMA}(\epsilon _1,\epsilon ;s),$$
(60)
and a part $`\delta v(\epsilon ,s)`$, that contains the deviations $`\stackrel{~}{F}`$. Self-consistency requires:

$$\delta v(\epsilon ,s)=0.$$
(61)
$`v(\epsilon )`$, as defined by Eq.(60), is in accordance with the definition (38), taking into account Eq.(B3) and the inequality $`\omega f^{}(\epsilon )/f(\epsilon )1`$ (Eq. (40)), owing to which contributions proportional to this parameter are negligible.

We now focus on the self-consistency equation (61). While Eq.(60) contains only the effective medium contribution to the diffusion propagator, $`\delta v`$ is a functional of the effective medium $`f`$ and the deviation $`\stackrel{~}{F}`$. By construction, it is at least linear in $`\stackrel{~}{F}`$. If this equation could be solved exactly, an exact solution to the diffusion problem could be found within quasi-elastic accuracy. In practise this is not possible, and, therefore, we depend on further approximations. To simplify this equation, we take into account only the lowest order contributions to this equation with respect to $`\stackrel{~}{F}`$, i.e. we linearize $`\delta v`$ with respect to $`\stackrel{~}{F}`$ and require that the first order contribution vanishes. This approach is quite close to the usual CPA-philosophy, in which vanishing of the t-matrix is required in its lowest order approximation. Using this procedure we obtain the following self-consistency equation:

$$\frac{1}{2}\omega ^2𝒩(ϵ)f(ϵ,s)\stackrel{~}{W}(0|ϵ;s)=a\omega \frac{sb}{\stackrel{~}{W}(0|ϵ;s)𝒩(ϵ)},$$
(62)
where $`a`$ and $`b`$ are simply numbers. A detailed derivation of this equation is given in Appendix C.

In deriving the self-consistency equations we have imposed further restrictions on the effective transition probabilities, which determine the range of its applicability. These restrictions can be formulated most conviniently using the dimensionless critical hopping length $`\rho _c(\epsilon ,s)`$, related to $`f(\epsilon ,s)`$ by (see Appendix B):

$$f(\epsilon ,s)\nu =\mathrm{exp}\rho _c(\epsilon ,s).$$
(63)
In terms of $`\rho _c`$, the inequalities (40) and (41), used in the derivation, read (prime is derivative with respect to $`\epsilon `$):

$$\left|\omega \rho _c^{}(\epsilon ,s)\right|1,$$
(64)
$$\left|\rho _c(\epsilon ,0)\rho _c(\epsilon ,s)\right|\rho _c(\epsilon ,0).$$
(65)
In addition, when calculating the integrals,

$$\rho _c(\epsilon ,s)1,$$
(66)
was used.

A closed solution to the self-consistency equation can only be found in the limit $`s=0`$. There we obtain:

$$\rho _c\left(\epsilon ,s=0\right)=\frac{2\alpha }{\left[\omega 𝒩\left(\epsilon \right)\right]^{1/d}}\left[\frac{2da}{S\left(d\right)}\right]^{1/d},$$
(67)
where $`S(d)`$ is the solid angle in $`d`$ dimensions.

For $`s`$ satisfying (65), Eq.(62) can be cast into the form:

$$\left[\rho _c(\epsilon ,0)\rho _c(\epsilon ,s)\right]\mathrm{exp}\left[\rho _c(\epsilon ,0)\rho _c(\epsilon ,s)\right]=\frac{s}{\mathrm{\Omega }\left(\epsilon \right)},$$
(68)
where:

$$\mathrm{\Omega }\left(\epsilon \right)=\frac{2da}{b\omega \rho _c(\epsilon ,0)}v(\epsilon ,0).$$
(69)
According to Eq.(68), the critical hopping length decreases with increasing frequency. Note, that the structure of the equation (68) for the calculation of the dispersion of the critical hopping length, obtained here, agrees completely with that obtained for the critical hopping length in calculating the equilibrium conductivity .

The dispersion of the transport coeffcicients is determined completely by the dispersion of the critical hopping length. For small $`s`$ the frequency dependence preexponetial factors can be ignored, so that from Eq.(68) explicit equations for the transport coefficients can be obtained. They are given by:

$$\frac{D(\epsilon ,s)}{D(\epsilon ,0)}\mathrm{ln}\frac{D(\epsilon ,s)}{D(\epsilon ,0)}=\frac{s}{\mathrm{\Omega }\left(\epsilon \right)},$$
(70)
$$\frac{v(\epsilon ,s)}{v(\epsilon ,0)}\mathrm{ln}\frac{v(\epsilon ,s)}{v(\epsilon ,0)}=\frac{s}{\mathrm{\Omega }\left(\epsilon \right)}.$$
(71)
The formal solution of these equations is given by the Lambert’s W-function $`𝒲(z)`$, defined by the equation $`z=𝒲(z)\mathrm{exp}𝒲(z)`$. Using the Lambert’s function we can write:

$$D(\epsilon ,s)=D(\epsilon ,0)\mathrm{exp}𝒲\left(s/\mathrm{\Omega }\left(\epsilon \right)\right),$$
(72)
$$v(\epsilon ,s)=v(\epsilon ,0)\mathrm{exp}𝒲\left(s/\mathrm{\Omega }\left(\epsilon \right)\right).$$
(73)
## VIII Constant density of states

A constant density of states, although of not much physical relevance, gives us the unique opportiunity to study pure dispersive energy transport. Here both $`\mathrm{\Omega }(\epsilon )`$ and $`v_0=v(\epsilon ,0)`$ are independent of energy. Consequently, Eq.(44) can readly be integrated. The integration yields:

$$P_L(\epsilon _0\epsilon ,s)=\frac{\mathrm{\Theta }\left(\epsilon _0\epsilon \right)}{v\left(s\right)}\mathrm{exp}\left[\frac{s\left(\epsilon _0\epsilon \right)}{v\left(s\right)}\right].$$
(74)
The time-dependence of this function can be obtained by the inverse Laplace-transformation. Using Eq.(71) to change the integration variable from $`s`$ to $`y=v/v_0`$, the inverse Laplace-transform of Eq.(74) may be written as:

$$P_L(\epsilon _0\epsilon ,t)=\frac{\mathrm{\Omega }}{v_0}\mathrm{\Theta }\left(\epsilon _0\epsilon \right)_C\frac{dy}{2\pi i}\frac{\mathrm{ln}y+1}{y}\mathrm{exp}\left[\mathrm{\Omega }ty\mathrm{ln}y\frac{\mathrm{\Omega }}{v_0}\left(\epsilon _0\epsilon \right)\mathrm{ln}y\right],$$
(75)
with the properly chosen integration contour $`C`$. At large enough $`t`$, that is for $`\mathrm{\Omega }t1`$, this expression, using the saddle-point method, gives simply the Gaussian packet:

$$P_L(\epsilon _0\epsilon ,t)\mathrm{\Theta }\left(\epsilon _0\epsilon \right)\frac{\mathrm{\Omega }}{v_0}\frac{1}{2\sqrt{\pi \mathrm{\Omega }t}}\mathrm{exp}\left[\frac{\left(\epsilon _0\epsilon v_0t\right)^2}{4v_0^2t/\mathrm{\Omega }}\right].$$
(76)
Now let us consider the time dependence of the energy relaxation. According to Eq.(51) the velocity of energy relaxation is:

$$\frac{dE(t)}{dt}=_C\frac{ds}{2\pi i}\frac{v\left(s\right)}{s}e^{st}.$$
(77)
Again, the time dependence of the energy relaxation speed can be calculated using asymptotics. Details of the calculations are presented in Appendix D, where it is shown that:

$$\frac{dE\left(t\right)}{dt}v_0\{\begin{array}{c}1+\sqrt{\frac{e}{2\pi }}\frac{1}{\left(\mathrm{\Omega }t\right)^{3/2}}\mathrm{exp}\left(\mathrm{\Omega }t/e\right),\mathrm{as}\mathrm{\Omega }t1,\hfill \\ \frac{\sqrt{2\pi }}{e\mathrm{\Omega }t}\left[\mathrm{ln}\left(e/(\mathrm{\Omega }t)\right)\right]^2,\mathrm{as}\mathrm{\Omega }t1.\hfill \end{array}$$
(78)
Note, that the problem of the energy relaxation in the case of constant density of states ($`v(\epsilon ,s)`$ is independent of $`\epsilon `$) is completely equivalent to one of the nonmarkovian charge transport in strong electric field $`E`$, when the diffusion (described by the second coordinate derivate) is totally neglected. One has only to replace $`v\left(s\right)u\left(s\right)E`$, $`u\left(s\right)`$ being the mobility. In the context of dispersive particle transport, the regimes $`\mathrm{\Omega }t1`$ and $`\mathrm{\Omega }t1`$ are the regimes of anomalous and normal “diffusion”, respectively. The main difference between both is in that while the sites are ususally distributed homogeneously in space, the density of states is usually an increasing function of energy.

## IX Exponential density of states

### A The saddle-point approximation and its break down for large times

The calculation of the time dependence of the energy distribution function for an arbitrary density of states on the basis of the Eqs. (36)-(38), (70) and (71) turns out to be a quite intricate problem. A tool that can be utilized in tackeling the problem is the saddle point approximation. How to apply the saddle point approximation for the calculation of the quantities of interest for an arbitrary density of states is shown in Appendix E. Below we focus on the exponential density of states, which is relevant, e.g., for amorphous Si:H.

We assume that the density of states is given by:

$$𝒩(\epsilon )=𝒩_0\mathrm{exp}\left(3\frac{\epsilon }{\mathrm{\Delta }}\right).$$
(79)
To simplify the notations we use the abbreviations:

$$\overline{\omega }=\frac{b\omega }{2da},$$
(80)
$$\overline{\nu }=\frac{2da^2}{b}\nu .$$
(81)
Then, the Eqs.(B3), (67) and (69) can be cast into the form:

$$v(\epsilon ,0)v\left(\epsilon \right)=\overline{\omega }\overline{\nu }\mathrm{exp}\left[\rho \left(\epsilon \right)\right],$$
(82)
$$\mathrm{\Omega }\left(\epsilon \right)=\frac{\overline{\nu }}{\rho \left(\epsilon \right)}\mathrm{exp}\left[\rho \left(\epsilon \right)\right],$$
(83)
$$\rho \left(\epsilon \right)=\rho _c(\epsilon ,0)=A\mathrm{exp}\left(\frac{\epsilon }{\mathrm{\Delta }}\right),$$
(84)
where:

$$A=\frac{2\alpha }{\left[\omega 𝒩\left(0\right)\right]^{1/d}}\left[\frac{2da}{S\left(d\right)}\right]^{1/d}.$$
(85)
Using these equations it follows from the formulas derived in Appendix E that the time dependence of the mean squared deviation and the time dependence of the mean energy are given by

$$\sigma ^2(\epsilon _m,\epsilon _0)\sigma ^2(t,\epsilon _0)\overline{\omega }\mathrm{\Delta }\left[1\frac{t_0^2}{\left(t+t_0\right)^2}\right],$$
(86)
$$\epsilon _m\left(t\right)\mathrm{\Delta }\mathrm{ln}\left[\frac{1}{A}\mathrm{ln}\left(\frac{\overline{\omega }}{\mathrm{\Delta }}\overline{\nu }\left(t+t_0\right)\right)\right],$$
(87)
where

$$t_0\left(\epsilon _0\right)=\frac{\mathrm{\Delta }}{\overline{\nu }\rho \left(\epsilon _0\right)v\left(\epsilon _0\right)}.$$
(88)
Note, that, according to the saddle-point approximation, the distribution is Gaussian. Furthermore, the dispersion is constant for $`tt_0`$. Consequently, for time scales in line with the applicability of the saddle-point approximation, the motion of the energy packet is “soliton-like”, that is the packet moves without distortion.

The applicability condition for the saddle-point (E13) requires:

$$\frac{\overline{\omega }}{\mathrm{\Delta }}\rho ^2(\epsilon _m)1.$$
(89)
Because of $`\rho (\epsilon _m)`$ grows with $`t`$, the saddle-point method, and, consequently, all the results of this chapter, becomes invalid, at least at sufficiently large $`t`$. The physical meaning of the condition (89) is the following: at $`tt_0`$ the width of the electron’s energies distribution is $`\delta \epsilon =\sqrt{\overline{\omega }\mathrm{\Delta }}`$. Because of $`\rho \mathrm{exp}\left(\epsilon /\mathrm{\Delta }\right)`$, the variation of $`\rho `$ is: $`\delta \rho /\rho =\sqrt{\overline{\omega }/\mathrm{\Delta }}`$, and, because $`v\mathrm{exp}\left(\rho \right)`$, we have: $`\delta v/v=\delta \rho =\sqrt{\overline{\omega }\rho ^2/\mathrm{\Delta }}`$. Thus, the condition (89) is just the one of the variation of electron’s “velocity” across the distribution is smaller, then the velocity itself.

### B Form of the distribution for large times

Even if the initial condition are such, that the condition (89) is fulfilled, and a gaussian distribution is formed, at some moment of time $`\rho _m\rho (\epsilon _m)`$ becomes of the order of $`\sqrt{\mathrm{\Delta }/\overline{\omega }}`$, which definitely should results in some deviation of the distribution function from its symmetric, gaussian form. For the exponential density of states, distribution function in the Laplace representation, Eq.(44), using Eqs.(73,82-84), may be written as:

$$P_L(\epsilon _0,\epsilon ;s)=\frac{1}{s}\rho \frac{}{\rho }\mathrm{exp}\left[\frac{1}{\overline{\omega }}\underset{\rho _0}{\overset{\rho }{}}𝑑x\frac{𝒲\left(sxe^x\right)}{x^2}\right]=\frac{𝒲\left(s\rho e^\rho \right)}{\overline{\omega }\rho s}\mathrm{exp}\left[\frac{1}{\overline{\omega }}\underset{\rho _0}{\overset{\rho }{}}𝑑x\frac{𝒲\left(sxe^x\right)}{x^2}\right],$$
(90)
where it was set $`\overline{\nu }=1`$ and $`\mathrm{\Delta }=1`$ by appropriate choice of time and energy units. Under the same conditions, which were used in the derivation of the above formula, the following approximation is valid (see Appendix F):

$$\underset{\rho _0}{\overset{\rho }{}}𝑑x\frac{𝒲\left(sxe^x\right)}{x^2}G(s,\rho )G(s,\rho _0);G(s,\rho )=\frac{1}{\rho ^2}𝒲\left(s\rho e^\rho \right)\left[1+\frac{1}{2}𝒲\left(s\rho e^\rho \right)\right]+O(\frac{𝒲}{\rho ^3},\frac{𝒲^3}{\rho ^3}).$$
(91)
When calculating distribution function at sufficiently large times, the characteristic values of $`s`$, giving the main contribution into inverse Laplace integral become so small, that one can suppose $`\left|s\right|\rho _0\mathrm{exp}(\rho _0)1`$, and: $`G(s,\rho _0)s\rho _0^1\mathrm{exp}\left(\rho _0\right)`$. Under this condition the distribution function approach initial conditions independent shape: $`P(\epsilon _0,\epsilon ;t)\varphi (\epsilon ,t+t_0(\epsilon _0))`$, $`t_0(\epsilon _0)=\left(\overline{\omega }\rho _0\right)^1\mathrm{exp}\rho _0`$,

$$\varphi (\epsilon ,t)=\underset{i\mathrm{}}{\overset{+i\mathrm{}}{}}\frac{ds}{2\pi i}\varphi _L(\epsilon ,s)e^{st}=\underset{i\mathrm{}}{\overset{+i\mathrm{}}{}}\frac{ds}{2\pi is}\rho \frac{}{\rho }\mathrm{exp}\left\{st\frac{1}{\overline{\omega }\rho ^2}𝒲\left(s\rho e^\rho \right)\left[1+\frac{1}{2}𝒲\left(s\rho e^\rho \right)\right]\right\}.$$
(92)
Let us consider, first, momenta of the distribution at the large times. In Laplace-representation the momenta are defined by the equation:

$$\chi _{nL}\left(s\right)=\underset{0}{\overset{\mathrm{}}{}}𝑑te^{st}\epsilon ^n\left(t\right)=\underset{\rho _0}{\overset{\mathrm{}}{}}\frac{d\rho }{\rho }\left(\mathrm{ln}\rho \right)^nP_L(\epsilon _0,\epsilon ,s)e^{st_0}\frac{d\rho }{\rho }\left(\mathrm{ln}\rho \right)^n\varphi _L(\epsilon ,s).$$
(93)
We shall omit further on irrelevant time shift multiple $`e^{st_0}`$. It turns out, that for $`\overline{\omega }\rho 1`$, the integrand of Eq.(93) has a maximum, located at:

$`\rho _m\left(1+\sqrt{\overline{\omega }}\right)\mathrm{ln}{\displaystyle \frac{\sqrt{\overline{\omega }}}{s}}.`$

However, if one would try to use the saddle-point method in evaluation of integrals like (90), the results would be incorrect at $`\sqrt{\overline{\omega }}\rho _m1`$, since the saddle point criterion is not fulfilled. However, at $`\overline{\omega }\rho ^21`$, momenta can be calculated as an expansion on the powers of small paramenter $`1/\left(\sqrt{\overline{\omega }}\rho \right)`$. Details of the calculations are presented in Appendix F.

As a result, we have for the first momentum, or the mean energy:

$`\chi _1\left(t\right)=\epsilon \left(t\right)\mathrm{\Delta }\mathrm{ln}{\displaystyle \frac{1}{A}}\mathrm{ln}\left(\stackrel{~}{b}\sqrt{\overline{\omega }}t\right)`$ (94)
$`\sqrt{{\displaystyle \frac{\pi \overline{\omega }}{2}}}+\sqrt{{\displaystyle \frac{\pi }{8\overline{\omega }}}}{\displaystyle \frac{1}{\mathrm{ln}^2\left(\stackrel{~}{b}\sqrt{\overline{\omega }}t\right)}}+\mathrm{},`$ (95)
where $`\stackrel{~}{b}=\sqrt{2}e^{1+\gamma /2}`$, and $`\gamma `$ is Eulers’s constant. For the distribution’s dispersion $`\sigma ^2(t)=\mu _2(t)=\chi _2(t)\chi _1^2(t)`$ we obtain:

$$\sigma ^2\left(t\right)=\mu _2\left(t\right)=\left(2\frac{\pi }{2}\right)\overline{\omega }+\left(1+\mathrm{ln}2\right)\frac{\sqrt{2\pi \overline{\omega }}}{\mathrm{ln}\left(\sqrt{\overline{\omega }}t\right)}+\mathrm{},$$
(96)
and for the third central momentum $`\mu _3(t)=\left(\epsilon \epsilon \left(t\right)\right)^3=\chi _3(t)3\chi _1\left(t\right)\chi _2\left(t\right)+2\chi _1^3\left(t\right)`$ we have:

$$\mu _3\left(t\right)=\left(\pi 3\right)\sqrt{\frac{\pi }{2}}\overline{\omega }^{3/2}+3\left[\pi \left(1\mathrm{ln}2\right)1\right]\frac{\overline{\omega }}{\mathrm{ln}\left(\sqrt{\overline{\omega }}t\right)},$$
(97)
Figures 5, 6, and 7 show mean energy, $`\chi _1`$, mean square deviation of energy $`\sigma `$, and dimensionless coefficient of asymmetry:

$$\mathrm{Asymmetry}=\frac{\mu _3\left(t\right)}{\sigma ^3\left(t\right)},$$
(98)
respectively, plotted versus $`\mathrm{ln}\mathrm{ln}t`$. There values, obtained numerically, using the distribution function (92), are compared with their corresponding asymptotic forms, following from Eqs.(95)-(97).

Moving along the same line, the whole distribution function may be reconstructed at large times, if to calculate all its momenta. However, it may be obtained in more simple way. Namely, when $`\overline{\omega }\rho ^21`$, one may replace $`s`$ within arguments of $`𝒲`$-functions in Eq.(92) with some $`s_0=c/t`$, $`c1`$ (a kind of argumentation may be found in Appendix G). Then the integral can easily be calculated to yield:

$`\varphi (\epsilon ,t)=`$ (99)
$``$ $`\rho {\displaystyle \frac{}{\rho }}\mathrm{exp}\left\{{\displaystyle \frac{1}{\overline{\omega }\rho ^2}}𝒲\left({\displaystyle \frac{c}{t}}\rho e^\rho \right)\left[1+{\displaystyle \frac{1}{2}}𝒲\left({\displaystyle \frac{c}{t}}\rho e^\rho \right)\right]\right\}.`$ (100)
One can easily show (see Appendix G), that asymptotic expressions for momenta (95-97) may be entirely reproduced from Eq.(100), if to set $`c=\mathrm{exp}(\gamma )`$ ($`\gamma `$ is the Euler’s constant). Figures 8-10 show some plots of the distribution function (92), compared with its asymptotic form (100). This latter can be further symplified, if to introduce $`\rho _t`$:

$$\frac{t}{c}=\rho _te^{\rho _t};\rho _t=𝒲\left(\frac{t}{c}\right).$$
(101)
As we shall see, the body of the distribution at a given time corresponds to values of $`\rho >\rho _t`$, $`\rho \rho _t\rho _t`$. Then, one can use approximation: $`𝒲\left(\left(c/t\right)\rho e^\rho \right)\rho \rho _t1`$. Then, we have in the exponent $`\left(1\rho _t/\rho \right)^2/2\overline{\omega }`$, which can be replaced by $`\left(\epsilon \epsilon _t\right)^2/2\overline{\omega }`$, where $`\epsilon _t\mathrm{ln}\mathrm{ln}t`$. As a result, at large time and small $`\overline{\omega }`$ the distribution function turns into:

$$\varphi (\epsilon ,t)=\frac{}{\epsilon }\{\begin{array}{c}\mathrm{exp}\left[\frac{\left(\epsilon \epsilon _t\right)^2}{2\overline{\omega }}\right],\epsilon <\epsilon _t\\ 1,\epsilon >\epsilon _t\end{array};\epsilon _t=\mathrm{ln}\mathrm{ln}t.$$
(102)
### C Initial stage of the evolution

At $`t=0`$ it was supposed $`P(\epsilon ,\epsilon _0;t=0)=\delta (\epsilon \epsilon _0)`$. To understand, what happens at the intitial stage of the distribution’s evolution, let us consider the behaviour of $`P`$ at $`\epsilon \epsilon _0`$. One can write from Eqs.(E1,73,82-84):

$`P_0(\epsilon _0,t)`$ $``$ $`P(\epsilon _0,\epsilon _0;t)={\displaystyle \frac{1}{\overline{\omega }\rho _0}}{\displaystyle \underset{i\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{ds}{2\pi is}}𝒲\left({\displaystyle \frac{s}{\overline{\nu }}}\rho _0e^{\rho _0}\right)e^{st}`$ (103)
$`=`$ $`{\displaystyle \frac{1}{\overline{\omega }\rho _0}}{\displaystyle _C}{\displaystyle \frac{dz}{2\pi i}}\left(1+z\right)\mathrm{exp}\left(\tau ze^z\right)`$ (104)
$`=`$ $`{\displaystyle \frac{1}{\overline{\omega }\rho _0\tau }}{\displaystyle _C}{\displaystyle \frac{dz}{2\pi i}}\mathrm{exp}\left(\tau ze^zz\right);`$ (106)
$`\tau ={\displaystyle \frac{\overline{\nu }t}{\rho _0}}e^{\rho _0},\rho _0=A\mathrm{exp}\left({\displaystyle \frac{\epsilon _0}{\mathrm{\Delta }}}\right),`$
where the new integration variable $`z=𝒲\left(\left(s/\overline{\nu }\right)\rho _0e^{\rho _0}\right)`$ was introduced instead of $`s`$, and the integration by parts was performed. This integral may be evaluated using the saddle-point method. The saddle-point equation, $`f^{}(z)=\tau (z+1)e^z1=0`$, gives us saddle point value $`z_c=𝒲\left(e/\tau \right)1`$, that is $`z_c\mathrm{ln}(1/\tau )\mathrm{ln}\mathrm{ln}(1/\tau )`$ at $`\tau 1`$, and $`z_ce/\tau 1`$ as $`\tau 1`$. To ensure the correctness of the saddle point approximation, one should require the parameter$`\left|f^{\prime \prime \prime }(z_c)\right|^2/\left|f^{\prime \prime }(z_c)\right|^3=\left|(3+z_c)^2/(2+z_c)^3\right|`$ to be small. While this is true at $`\tau 1`$, this parameter appears to be $`1`$ at $`\tau 1`$. Therefore, the result for $`\tau 1`$ is correct up to the multiple of the order of 1 only. We have:

$$P_0(\rho ,t)\{\begin{array}{c}\frac{1}{\sqrt{2\pi }\overline{\omega }\rho }\left(\frac{e}{\tau }\right)^{3/2}\mathrm{exp}\left(\frac{\tau }{e}\right)\text{ as }\tau 1,\\ \frac{e}{\sqrt{2\pi }\overline{\omega }\rho }\left(\mathrm{ln}\frac{e}{\tau }\mathrm{ln}\mathrm{ln}\frac{e}{\tau }\right)\text{as }\tau 1.\end{array}$$
(107)
In a similar fasion one can calculate the energy derivative at $`\epsilon =\epsilon _0`$:

$`P_1(\rho _0,t){\displaystyle \frac{P}{\epsilon }}|_{\epsilon =\epsilon _0}={\displaystyle \frac{\rho }{\mathrm{\Delta }}}{\displaystyle \frac{}{\rho }}P(\rho _0,\rho ;t)|_{\rho =\rho _0}`$ (108)
$`=`$ $`{\displaystyle \frac{1}{\overline{\omega }\mathrm{\Delta }\rho _0}}{\displaystyle \underset{i\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{ds}{2\pi is}}𝒲\left({\displaystyle \frac{\mathrm{\Delta }𝒲}{\overline{\omega }\rho _0}}+{\displaystyle \frac{𝒲\rho _0}{1+𝒲}}\right)e^{st}`$ (109)
$``$ $`{\displaystyle \frac{1}{\left(\overline{\omega }\rho _0\right)^2}}{\displaystyle _C}{\displaystyle \frac{dz}{2\pi i}}\left[z\left(1+z\right){\displaystyle \frac{\overline{\omega }}{\mathrm{\Delta }}}\rho _0^2\right]\mathrm{exp}\left(\tau ze^z\right),`$ (110)
where the term $`(\overline{\omega }/\mathrm{\Delta })\rho _0z`$ within the square brackets was neglected due to assumption $`\left(\overline{\omega }/\mathrm{\Delta }\right)\rho 1`$. At $`\tau 1`$ this integral may be readily evaluated in the saddle point approximation, which yields:

$$P_1(\rho ,t)\frac{1}{\overline{\omega }\mathrm{\Delta }}\sqrt{\frac{e}{2\pi \tau }}\mathrm{exp}\left(\frac{\tau }{e}\right).$$
(111)
When $`\tau 1`$, one have to perform integration by parts first, which gives:

$`P_1(\rho ,t)={\displaystyle \frac{1}{\overline{\omega }^2\rho ^2\tau }}`$
$`{\displaystyle _C}{\displaystyle \frac{dz}{2\pi i}}\left[z1{\displaystyle \frac{\overline{\omega }\rho ^2}{\mathrm{\Delta }\left(1+z\right)}}{\displaystyle \frac{\overline{\omega }\rho ^2}{\mathrm{\Delta }\left(1+z\right)^2}}\right]\mathrm{exp}\left(\tau ze^zz\right).`$
Then one have in the saddle point approximation:

$`P_1(\rho ,t)`$ $``$ $`{\displaystyle \frac{e}{\sqrt{2\pi }\overline{\omega }^2\rho ^2}}\left(z_c^21\overline{\omega }\rho ^2\right),`$ (112)
$`z_c`$ $``$ $`\mathrm{ln}{\displaystyle \frac{1}{\tau }}\mathrm{ln}\mathrm{ln}{\displaystyle \frac{1}{\tau }}.`$ (113)
At sufficiently small $`t`$ the above expression is positive, which means that the distribution is “sticked” near $`\epsilon =\epsilon _0`$, that is, it monotoneously decreases as $`\epsilon _0\epsilon >0`$ increases. At some $`\tau =\tau _0`$, $`P_1`$ changes its sign. If $`\left(\overline{\omega }/\mathrm{\Delta }\right)\rho _0^21`$, this corresponds to $`z_{c0}\mathrm{ln}\frac{1}{\tau _0}\sqrt{\overline{\omega }/\mathrm{\Delta }}\rho _0`$, otherwise $`\tau _01`$. This correspons just to the moment of time, when this distribution separates from the intial point, — its maximum is at $`\epsilon <\epsilon _0`$. For $`\tau \tau _0`$ the width of the distribution on the energy scale may be estimated as:

$$\mathrm{\Delta }\epsilon (t)\frac{P_0(\rho ,t)}{P_1(\rho ,t)}\frac{\overline{\omega }\rho }{z_c}\frac{\overline{\omega }\rho }{\mathrm{ln}(1/\tau )}.$$
(114)
This formula is valid if $`\mathrm{ln}(1/\tau )\rho `$, up to the times corresponding either $`\tau 1`$ if $`\left(\overline{\omega }/\mathrm{\Delta }\right)\rho _0^21`$, or $`\mathrm{ln}\left(1/\tau \right)\sqrt{\overline{\omega }/\mathrm{\Delta }}\rho _0`$ otherwise. In the former case the distribution width becomes of the order of $`\overline{\omega }\rho `$ just before it separates from the initial point. After this, the distribution, as it was shown in previous subsection, widens till $`\sqrt{\overline{\omega }\mathrm{\Delta }}`$, and the gaussian packet moves downward, untill its center reaches the value, corresponding to $`\rho \sqrt{\mathrm{\Delta }/\overline{\omega }}`$. In the latter case the packet’s width is $`\sqrt{\overline{\omega }\mathrm{\Delta }}\overline{\omega }\rho `$ at the very moment of its “unsticking”. In both cases, to consider subsequent evolution, one has to use some other approximation instead of saddle point one.

## X Conclusions

In the paper we have presented an effective method which permitts the investigation of relaxation phenomena of localized charge carriers far from equillibrium due to phonon-assisted hopping at zero temperature. From the point of view of the formalism the main equations are the equations for the calculation of the diffusion propagator (36), (37) and (38), and the equations (69), (70) and (71), which determine the dispersion of the transport coefficients. These equations show, that both particle transport and energy transport are dispersive. The equations (70) and (71) lead to a strong dependence of the diffusion constant $`D(ϵ;s)`$ and the velocity of energy relaxation $`v(ϵ;s)`$ on frequency $`s`$ already for low frequencies. To our knowledge, these equations have not been dervied in the literature so far, for systems far from equillibrium. In fact, in the literature mainly frequency independent transport coefficients can be found (see e.g.  and references therein). The strong dependence of the transport coefficients on $`s`$ results in a non-Markovian equation for the calculation of the diffusion propagator (36). This discriminates our equations from the Markovian integral equations used in the Refs. -. The latter equations can, in principle, be obtained from our leading equation (11) by neglecting statistical correlation, so as to average every factor independently.

Using our effective-medium method, we have investigated the relaxation of charge carriers in band tails at $`T=0`$. According to our results, energy relaxation is connected with dispersional transport. Even if we have no real diffusion in the system, since the temperature is zero, we have some spreading of the energy distribution with time. For a constant density of states the situation is completely equivalent to the one in the problem of the electron’s motion in a disordered system, subjected to an electric field . The main difference between the above mentioned problem and the energy relaxation is in that, while in the former case the sites are distributed homogeneously in space, in the energy relaxation problem the density of states is a decaying function of energy for most physical systems. In this situation we arrive at the picture that the particles of a package with lower energies move slower than particles with higher energy. This leads to the opposite tendency: at first, there is a slowing of the package spreading. The time dependence of the dispersion, being linear at the first stage of the evolution, is later slowing down. Later on, different possibilities exist, dependend on the particular type of the energy dependence of the density of states.

We performed a detailed investigation for the exponential density of states for two time regimes. If the variation of the velocity $`v`$ in the energy space across the distribution is smaller than the velocity itself, that is if $`\omega \rho ^2/\mathrm{\Delta }1`$, we arrive at the situation of a package, moving steadily down along the energy axis without any deformation. Dispersion becomes time independent (see Eq.(86)). However, since the particles are sinking down this condition becomes violated when time goes by. The parameter $`\omega \rho ^2/\mathrm{\Delta }`$, being small at the first stage of the evolution, is getting large. When this parameter is larger than one, the steady motion condition becomes violated again. The package, being of gaussian form before, undergoes some restructuring to another, non-gaussian stable form, with its width lower than before by some numeric factor of the order of 1 (Eq.(96)). In our pictures, Figs 7-10, it is clearly seen, that the package becomes non-gaussian. This result remains valid until the very moment the quasi-elasticity conditions break down ($`\omega \rho /\mathrm{\Delta }1`$).

For the exponential density of states we have found the packet to move as $`\mathrm{ln}\mathrm{ln}t`$, i.e., its motion is strongly slowing down with time, the packet becomes, roughly speaking, almost stopped. This type of behaviour may be called “glassy”, because of the overall time scale for the packet evolution (governed by the exponential function of the large parameter $`\mathrm{\Delta }/\omega `$), becomes huge. For example, when $`\omega /\mathrm{\Delta }=10^4`$, this time scale approaches $`10^{1000}`$ and more (see Fig. 10)!

The main simplification used in our paper is the quasi-elastic approximation. This approximation relies on the smallness of the upper bound of the energy transfer to the phonon system by one hop. For localized electrons, this upper bound can be much smaller as the Debye energy of the host material, since not all phonons can interact with localized electrons equally well. For localized electrons the electron phonon coupling constant approaches zero, for phonons with wave-vector $`q>2\alpha `$. Thus high energetic phonons are less effective. Only phonons with energies $`\omega <\omega _D2\alpha a`$, where $`a`$ is the lattice constant of the host material, are effective. Furthermore, in disordered systems the high energetic acustical phonons are localized, and thus need not contribute to transport by necessity. Nevertheless, the question whether the quasi-elastic approximation is applicable depends much on the material of interest. If, however, we compare our result for the time dependence of the mean energy (87) in an exponential density of states with those existing in the literature , we also find agreement for $`\omega /\mathrm{\Delta }1`$, which indicates that our results are, at least qualitatively, of wider validity. Unfortunately, to our knowledge, in the literature there are no further results on the width of the energy distribution available to compare with.

From our point of view the main question remained open is how the results, obtained for the exponential density of states, may be generalized for other types of energy dependencies. One can imagine, e.g., that for densities of states decaying with decreasing energy slower than exponential, the dispersion is growing, but under some sublinear law. Also, in the opposite case, in a density of states decaying faster than exponential, we would possibly have a dispersion tending to zero for large times.

It should, however, also be mentioned that it is not completely clear under which situation zero temperature results can be applied to systems at finite temperature. In our approach we assumed that $`T=0`$. However, in a real system, when the carriers are sinking down, the criterion for the temperature to be treated as zero, is violated at every finite value of $`T`$ at some moment of time, since the contribution of hops to sites with higher energy values becomes more and more comparable with one of downward hops. Therefore, the consideration of temperature becomes vital. Work in this direction is in progress.

## ACKNOWLEDGMENTS

On of the authors (O.B.) would like to thank Prof. P. Thomas for drawing his attention to the photoluminscense problem and stimulating discussions.

## A On the configuration average

Here we derive the set of integral equations (17)-(22). The same method has already been used in Ref..

According to Eq.(10), the configuration average of the electron density $`n(\rho ,s)`$ is given by:

$$n(\rho ,s)=𝑑\rho _1p_0\left(\rho _1\right)P(\rho _1,\rho ),$$
(A1)
where

$$P(\rho ^{},\rho )=\eta \left(\rho ^{}\right)\mathrm{\Phi }(\rho ^{},\rho )$$
(A2)
is the diffusion function. According to Eq.(11), the diffusion function is given by:

$$P=\frac{1}{s}\underset{n=0}{\overset{\mathrm{}}{}}\eta \left(\frac{V}{s}\right)^n.$$
(A3)
Diagrammatically this series can be represented as depicted in Fig.1. In this picture every full dot represents a potential $`V`$, which depends on disorder via the structur factor $`\eta `$. The expansion starts with a single structural factor $`\eta `$, depicted by the empty dot. Note that all disorder is comprised into the structure factor $`\eta `$.

In order to calculate the configuration average, we average the series (A3) term by term. The resulting expansion can be depicted in the usual form, as in impurity scattering problem (see Fig.2). One can see from the picture, that the set of all diagrams can be decomposed into two subsets, the set $`S`$ containing all diagrams connected with the empty point and the set $`F`$ containing all other diagrams. Owing to this decomposition configuration averaged diffusion function can be written as:

$$P(\rho ^{},\rho )=𝑑\rho _1S(\rho ^{},\rho _1)F(\rho _1,\rho ),$$
(A4)
from which Eq.(17) follows. Note that $`F`$ is nothing but $`\mathrm{\Phi }`$.

Again the class of diagrams contributing to $`F`$ can be decomposed into reducible and irreducible diagrams. All irreducible diagrams can be comprised into a function $`\mathrm{\Pi }`$. Doing so, we obtain Eq.(18). Note, that, although the introduction of the irreducible part $`\mathrm{\Pi }`$ is parallel to the introduction of the self-energy it has to be stressed, that its physical interpretation is completely different. In the present context $`\mathrm{\Pi }(\rho ^{},\rho )`$ is the true transition probability from $`\rho ^{}`$ to $`\rho `$. It is $`\mathrm{\Pi }(\rho ^{},\rho )`$ that has to be calculated and compared with the experimental situation and not the bar transition probability, as suggested in many papers (see e.g. ). Diagrams, contributing to $`\mathrm{\Pi }`$ are depicted in Fig.3. The set of diagrams depicted in Fig. 3 can be generated by means of the propagator $`\mathrm{\Pi }_{\stackrel{~}{\rho }}`$, defined by the equation (21), using Eq.(19).

The propagator $`\mathrm{\Pi }_{\stackrel{~}{\rho }}`$ can also be used to generate the irreducible block $`S`$. Diagrams contributing to $`S`$ are depicted in Fig.4. From the picture it follows, that $`S`$ is given by Eq.(22).

## B Transport coefficients

According to Eq.(38), $`v(\epsilon ,s)`$ is given by:

$$v(\epsilon ,s)=\frac{1}{2}\omega ^2𝒩\left(\epsilon \right)\stackrel{~}{W}\left(q=0|\epsilon \right)=\frac{1}{2}\omega ^2𝒩\left(\epsilon \right)\nu \mathrm{exp}\left[\rho _c(\epsilon ,s)\right]𝑑𝐑\frac{1}{1+\mathrm{exp}\left[2\alpha R\rho _c(\epsilon ,s)\right]},$$
(B1)
where:

$$f(\epsilon ,s)\nu =\mathrm{exp}\rho _c(\epsilon ,s).$$
(B2)
We assume, that $`\rho _c`$ is large, so the integrals can be calculated in the limit $`\rho _c\mathrm{}`$. This assumption is justified a posteriory by explicit calculation of $`\rho _c`$. In limit $`\rho _c\mathrm{}`$, the integrand reduces to a step function, so that the integrals can be calculated easily. Doing so, we find:

$$v(\epsilon ,s)=\frac{1}{2}\frac{S\left(d\right)}{d}\omega ^2𝒩\left(\epsilon \right)\left[\frac{\rho _c(\epsilon ,s)}{2\alpha }\right]^d\nu \mathrm{exp}\left[\rho _c(\epsilon ,s)\right],$$
(B3)
where $`S(d)`$ is the solid angle and $`d`$ is the spatial dimension. From the calculation and the result it is to be seen, that $`\rho _c(\epsilon ,s)`$ is to be identified with the dimensionless critical hopping length at energy $`\epsilon `$.

The same procedure is applied to the calculation of $`D(\epsilon ,s)`$, defined by Eq.(37). Using the same approximations, we obtain:

$$D(\epsilon ,s)=\frac{1}{2}\frac{S\left(d\right)}{d\left(d+2\right)}\omega 𝒩\left(\epsilon \right)\left[\frac{\rho _c(\epsilon ,s)}{2\alpha }\right]^{d+2}\nu \mathrm{exp}\left[\rho _c(\epsilon ,s)\right].$$
(B4)
Thus $`D(\epsilon ,s)`$ and $`v(\epsilon ,s)`$ differ only in preexponential factors.

We anticipate that for $`s=0`$ we obtain:

$$\rho _c(\epsilon ,0)=\frac{2\alpha }{\left[\omega 𝒩\left(\epsilon \right)\right]^{1/d}}\left[\frac{2da}{S\left(d\right)}\right]^{1/d},$$
(B5)
where $`a`$ is a number and $`S(d)`$ is the solid angle. A different form for the diffusion coefficient can be found in Ref.. The expression given in this paper differs from our expression essentially only in that $`\omega `$ is replaced by the tailing parameter. In looking on this difference one should, however, take into account, that the model used in that paper is different. Furthermore, in Ref. the transport coefficient are independent of $`s`$. It is therefore not completely clear to which range of time they apply. Clearly, they can not determine the evolution for all times.

## C Derivation of self-consistency equation

Here we derive the self-consistency equation (62). In order to work out explicitely the first order contribution to the self-consistency equation (61) we first need the first order correction to $`\mathrm{\Pi }`$. It is given by:

$$\mathrm{\Pi }^{(1)}(\rho ^{},\rho ;s)=𝑑\rho _1𝑑\rho _2𝑑\rho _3𝒩\left(\epsilon _3\right)\stackrel{~}{w}_{\rho _3}(\rho ^{},\rho _1;s)\stackrel{~}{F}(\rho _1,\rho _2)\stackrel{~}{w}_{\rho _3}(\rho _2,\rho ;s).$$
(C1)
If we insert the expression for $`\stackrel{~}{w}_{\rho _3}`$, we obtain:

$`\mathrm{\Pi }^{(1)}(\rho ,\rho ;s)={\displaystyle }d\rho _1`$ (C2)
$`\{𝒩\left(\epsilon \right)\stackrel{~}{W}(\rho ,\rho ;s)[\stackrel{~}{F}(\rho ,\rho _1)\stackrel{~}{F}(\rho ,\rho _1)]\stackrel{~}{W}(\rho _1,\rho ;s)𝒩\left(\epsilon _1\right)\stackrel{~}{W}(\rho ^{},\rho _1;s)[\stackrel{~}{F}(\rho _1,\rho )\stackrel{~}{F}(\rho ^{},\rho )]\stackrel{~}{W}(\rho ,\rho _1;s)\}.`$ (C3)
Thus, at this stage the self-consistency equation takes the form:

$`0`$ $`=`$ $`{\displaystyle }d𝐪d\epsilon d\epsilon _1(\epsilon \epsilon )\{𝒩\left(\epsilon \right)\stackrel{~}{W}(0|\epsilon ,\epsilon ;s)\stackrel{~}{F}(q|\epsilon ,\epsilon _1)\stackrel{~}{W}(q|\epsilon _1,\epsilon ;s)𝒩\left(\epsilon \right)\stackrel{~}{W}(q|\epsilon ^{},\epsilon ;s)\stackrel{~}{F}(q|\epsilon ^{},\epsilon _1)\stackrel{~}{W}(q|\epsilon _1,\epsilon ;s)`$ (C5)
$`𝒩\left(\epsilon _1\right)\stackrel{~}{W}(0|\epsilon ^{},\epsilon _1;s)\stackrel{~}{F}(q|\epsilon _1,\epsilon )\stackrel{~}{W}(q|\epsilon ,\epsilon _1;s)+𝒩\left(\epsilon _1\right)\stackrel{~}{W}(q|\epsilon ^{},\epsilon _1;s)\stackrel{~}{F}(q|\epsilon ^{},\epsilon )\stackrel{~}{W}(q|\epsilon ,\epsilon _1;s)\}.`$
In order to eliminate $`\stackrel{~}{F}`$, we have to replace $`\stackrel{~}{F}(q|\epsilon ^{},\epsilon )=F(q|\epsilon ^{},\epsilon )f(\epsilon ,s)\delta (\epsilon ^{}\epsilon )`$, where $`F`$ is the effective medium approximation of the diffusion propagator. Let us first work out the local contribution. Replacing $`\stackrel{~}{F}`$ by $`f\delta `$, the first and the third term of Eq.(C5) cancel each other. The 4th term is zero, taken into account that the delta function of the effective medium is multiplied by its argument. Thus, when $`\stackrel{~}{F}`$ is replaced by $`f(\epsilon ,s)\delta (\epsilon ^{}\epsilon )`$, only the 2nd term of Eq.(C5) survives. To simplify this term we take into account that

$$\stackrel{~}{W}(q|ϵ;s)=\stackrel{~}{W}(0|ϵ,s)\varphi (\frac{q\rho _c(ϵ,s)}{2\alpha }),$$
(C6)
where $`\varphi `$ is a dimensionless function. Consequently, we obtain

$``$ $`{\displaystyle 𝑑𝐪𝑑\epsilon ^{}\left(\epsilon ^{}\epsilon \right)𝒩\left(\epsilon \right)f(\epsilon ^{},s)\left[\stackrel{~}{W}\left(q|\epsilon ^{},\epsilon ;s\right)\right]^2}{\displaystyle 𝑑𝐪\left[\stackrel{~}{W}\left(𝐪|\epsilon ,s\right)\right]^2\frac{\omega ^2}{2}f(\epsilon ,s)𝒩\left(\epsilon \right)}`$ (C7)
$`=`$ $`{\displaystyle \frac{1}{2}}\omega ^2\left({\displaystyle \frac{2\alpha }{\rho _c(ϵ,s)}}\right)^d𝒩(ϵ)f(ϵ,s)\stackrel{~}{W}^2(0|ϵ,s){\displaystyle d^dx\varphi ^2(x)}.`$ (C8)
Here terms proportional to $`\omega f^1(\epsilon ,s)df(\epsilon ,s)/d\epsilon 1`$ have been neglected.

Let us now focus on the contribution of the regular part of the diffusion propagator to the self-consistency equation. Owing to the step functions in $`F`$ an $`\stackrel{~}{W}`$, the first and third term in Eq.(C5) are zero, when $`\stackrel{~}{F}`$ is replaced by $`F`$. Thus we are left with:

$$d𝐪d\epsilon d\epsilon _1(\epsilon \epsilon )\{𝒩\left(\epsilon \right)\stackrel{~}{W}(q|\epsilon ,\epsilon ;s)F(q|\epsilon ,\epsilon _1)\stackrel{~}{W}(q|\epsilon _1,\epsilon ;s)+𝒩\left(\epsilon _1\right)\stackrel{~}{W}(q|\epsilon ^{},\epsilon _1;s)F(q|\epsilon ^{},\epsilon )\stackrel{~}{W}(q|\epsilon ,\epsilon _1;s)\}.$$
(C9)
The range of integration in (C9) is determined by the step functions in $`F`$ and $`\stackrel{~}{W}(q|\epsilon ^{},\epsilon ;s)=\theta (\epsilon ^{}\epsilon )\theta (\omega \epsilon ^{}+\epsilon )\stackrel{~}{W}(q|\epsilon ^{};s)`$. Owing to these step functions, the energy integrations extend at most over intervals of length $`\omega `$, so that the quasi-elastic approximation can again be applied to the effective-transition probabilities and the density of states. The diffusion propagator enetering Eq.(C9), however, can not be dealt with in this way since for arbitrary $`q`$ the derivatives of the diffusion propagator are not small as compared to the diffusion propagator itself. Therefore, another procedure is needed. To simplify this expression further we consider Eq.(31). If terms small with respect to $`\omega N^{}/N`$ and $`\omega \stackrel{~}{W}^{}(q|ϵ,s))/\stackrel{~}{W}(q|ϵ,s)`$ are neglected the function

$$\mathrm{\Phi }(x|y^{},y)=\omega ^2F(x\frac{2\alpha }{\rho _c(\omega y,s)}|\omega y^{},\omega y)\stackrel{~}{W}(\omega y,s)𝒩(\omega y),$$
(C10)
can be introduced, that satisfies the equation

$`{\displaystyle \frac{s}{\stackrel{~}{W}(\omega y,s)𝒩(\omega y)\omega }}\mathrm{\Phi }(x|y^{},y)=\delta (y^{}y)+{\displaystyle \underset{0}{\overset{1}{}}}𝑑y_1[\mathrm{\Phi }(x|y^{},y_1+y)\varphi (x)\mathrm{\Phi }(x|y^{},y)\varphi (0)].`$ (C11)
Then the expression (C9) can be cast into the form

$$\omega \stackrel{~}{W}(\omega y)(\frac{2\alpha }{\rho _c(\omega y,s)})^dd^dx\varphi ^2(x)\left[\underset{0}{\overset{1}{}}𝑑y^{}\underset{0}{\overset{y^{}}{}}𝑑y_1y^{}\mathrm{\Phi }(x|y^{}+\frac{ϵ}{\omega },y_1+\frac{ϵ}{\omega })+\underset{0}{\overset{1}{}}\underset{0}{\overset{1y_1}{}}𝑑y^{}y^{}\mathrm{\Phi }(x|y^{}+\frac{ϵ}{\omega },\frac{ϵ}{\omega })\right]$$
(C12)
For $`s=0`$ Eq.(C11) does not contain any physical parameter. It only leads to the determination of the function

$`\mathrm{\Phi }_0(x|y^{},y)=\theta (y^{}y)\mathrm{\Phi }_0(x|y^{}y),`$

which satisfies the equation

$`0=\delta (y^{}y)+{\displaystyle \underset{0}{\overset{1}{}}}𝑑y_1[\mathrm{\Phi }_0(x|y^{},y_1+y)\varphi (x)\mathrm{\Phi }_0(x|y^{},y)\varphi (0)].`$ (C13)
Provided, we restrict our consideration to small frequencies in deriving the self-consistency equation we only need to take into account the linear contribution of the function $`\mathrm{\Phi }`$ with respect to $`s/(\stackrel{~}{W}(\omega y)𝒩(\omega y)\omega )`$. Then, using again the smallness of the variation of the effective transition probabilities and the density of states with respect to changes of energy over intervals of length $`\omega `$, $`\mathrm{\Phi }`$ can be approximated as

$$\mathrm{\Phi }(x|y^{},y)=\mathrm{\Phi }_0(x|y^{}y)\frac{s}{\stackrel{~}{W}(\omega y)𝒩(\omega y)\omega }\underset{y}{\overset{y^{}}{}}𝑑y_1\mathrm{\Phi }_0(x|yy_1)\mathrm{\Phi }_0(x|y_1y).$$
(C14)
Using this expression the self-consistency equation takes the form (62), where the coefficients $`a`$ and $`b`$ are given by

$$a=\frac{1}{2}\frac{d^dx\varphi ^2(x)\underset{0}{\overset{1}{}}(1y)^2\mathrm{\Phi }_0(x|y)}{d^dx\varphi ^2(x)},$$
(C15)
$$b=\frac{1}{2}\frac{d^dx\varphi ^2(x)_0^1𝑑y(1y)^2_0^y𝑑y_1\mathrm{\Phi }_0(x|yy_1)\mathrm{\Phi }_0(x|y_1)}{d^dx\varphi ^2(x)}.$$
(C16)
## D Calculation of the energy relaxation speed for constant DOS

To calculate the integral (77), we again use Eq.(71) to change the integration variable from $`s`$ to $`y=v/v_0`$. Furthermore, for convinience, we introduce the parameter $`\tau =\mathrm{\Omega }t`$. Doing so, we obtain:

$$\frac{dE\left(t\right)}{dt}=v_0F\left(\tau \right),$$
(D1)
where $`F(\tau )`$ is given by:

$$F\left(\tau \right)=_C\frac{dy}{2\pi i}\frac{1+\mathrm{ln}y}{\mathrm{ln}y}\mathrm{exp}\left(\tau y\mathrm{ln}y\right).$$
(D2)
After an integration by parts we obtain:

$$\frac{dF}{d\tau }=\frac{1}{\tau }_C\frac{dy}{2\pi i}\mathrm{exp}\left(\tau y\mathrm{ln}y\right)\frac{1}{\tau }S\left(\tau \right).$$
(D3)
To calculate $`F(\tau )`$, Eq. (D3) should be supplemented with the condition $`F(+\mathrm{})=1`$.

To evaluate $`S\left(\tau \right)`$, it is convenient to use the integration contour $`Im(y\mathrm{ln}y)=0`$, or, introducing polar coordinates $`y=r\mathrm{exp}\left(i\theta \right)`$:

$$r=\mathrm{exp}\left(\theta \mathrm{cot}\theta \right).$$
(D4)
Substituting Eq.(D4) into Eq.(D3), we have:

$$S\left(\tau \right)=_\pi ^\pi \frac{d\theta }{2\pi }V\left(\theta \right)e^{\tau V\left(\theta \right)},V\left(\theta \right)=\frac{\theta }{\mathrm{sin}\theta }e^{\theta \mathrm{cot}\theta }.$$
(D5)
The value of the integral may be estimated by the saddle-point method, looking for the maxima of the expression $`\mathrm{ln}V\left(\theta \right)\tau V\left(\theta \right)`$. The stationary point equation is:

$$V^{}\left(\theta \right)\left[\tau \frac{1}{V\left(\theta \right)}\right]=0.$$
(D6)
As $`\tau 1`$, the stationary point is $`\theta _s=0`$, and the asymptote of Eq.(D5) is:

$$S\left(\tau \right)\left(2\pi e\tau \right)^{1/2}\mathrm{exp}\left(\tau /e\right).$$
(D7)
On the other hand, if $`\tau 1`$, we have two stationary points $`\pm \theta _s`$, $`V\left(\pm \theta _s\right)=1/\tau `$, $`\theta _s=\pi \delta `$, and for $`\delta 1`$ the following equation may be obtained:

$`{\displaystyle \frac{\pi }{\delta }}\mathrm{exp}\left({\displaystyle \frac{\pi }{\delta }}1\right)={\displaystyle \frac{1}{\tau }},`$

the solution for which at small $`\tau `$ is: $`\pi /\delta =𝒲(e/\tau )\mathrm{ln}\left(e/\tau \right)\mathrm{ln}\mathrm{ln}(e/\tau ).`$The asymptotics of $`S\left(\tau \right)`$ turns out to be:

$$S\left(\tau \right)\frac{\sqrt{2\pi }}{e\tau }\left[\mathrm{ln}\left(e/\tau \right)\right]^2.$$
(D8)
## E The saddle-point approximation

For energy dependent density of states it turns out to be difficult to obtain explicite expressions for the time dependence of the energy distribution function. Here, according to Eq.(44), the time-dependence has to be calculated from the equation:

$$P(\epsilon ,\epsilon _0,t)=\underset{i\mathrm{}}{\overset{+i\mathrm{}}{}}\frac{ds}{2\pi iv(\epsilon ,s)}\mathrm{exp}\left[sts\underset{\epsilon }{\overset{\epsilon _0}{}}\frac{d\epsilon ^{}}{v(\epsilon ^{},s)}\right],$$
(E1)
where $`v(\epsilon ,s)`$ is given by Eq. (71) or (73). Let us try to integrate over $`s`$, using the saddle-point approximation. The saddle-point position $`s_0(\epsilon ,\epsilon _0,t)`$ may be obtained from the equation:

$$t=\underset{\epsilon }{\overset{\epsilon _0}{}}\frac{d\epsilon ^{}}{v\left(\epsilon ^{}\right)g(\epsilon ^{},s_0)}s_0\underset{\epsilon }{\overset{\epsilon _0}{}}\frac{d\epsilon ^{}}{v\left(\epsilon ^{}\right)g^2(\epsilon ^{},s_0)}\frac{g(\epsilon ,s_0)}{s_0},$$
(E2)
where $`g=v(\epsilon ,s)/v(\epsilon ,0)=\mathrm{exp}𝒲(s/\mathrm{\Omega }(\epsilon ))`$ was introduced. The expression for the diffusion propagator in the saddle-point approximation may be written as:

$$P_L(\epsilon ,\epsilon _0,t)=\frac{1}{\sqrt{4\pi D(\epsilon ,\epsilon _0,s_0)}}\frac{1}{v_0\left(\epsilon \right)g(\epsilon ,s_0)}\mathrm{exp}\left[s_0^2\underset{\epsilon }{\overset{\epsilon _0}{}}\frac{d\epsilon ^{}}{v_0\left(\epsilon ^{}\right)g^2(\epsilon ^{},s_0)}\frac{g(\epsilon ,s_0)}{s_0}\right],$$
(E3)
where:

$$D(\epsilon ,\epsilon _0,s_0)=\frac{^2}{s_0^2}\underset{\epsilon }{\overset{\epsilon _0}{}}\frac{d\epsilon ^{}\mathrm{\Omega }\left(\epsilon ^{}\right)\mathrm{ln}g(\epsilon ^{},s_0)}{2v_0\left(\epsilon ^{}\right)}>0.$$
(E4)
Assuming $`s_0/\mathrm{\Omega }1`$, we have: $`g(\epsilon ,s_0)1+s_0/\mathrm{\Omega }s_0^2/2\mathrm{\Omega }^2`$,

$$s_0=\frac{tT(\epsilon ,\epsilon _0)}{2D(\epsilon ,\epsilon _0)},D(\epsilon ,\epsilon _0)=\underset{\epsilon }{\overset{\epsilon _0}{}}\frac{d\epsilon ^{}}{v_0\left(\epsilon ^{}\right)\mathrm{\Omega }\left(\epsilon ^{}\right)},T(\epsilon ,\epsilon _0)=\underset{\epsilon }{\overset{\epsilon _0}{}}\frac{d\epsilon ^{}}{v_0\left(\epsilon ^{}\right)}$$
(E5)
and, finally, the diffusion propagator:

$$P_L(\epsilon ,\epsilon _0,t)\frac{1}{\sqrt{4\pi D(\epsilon ,\epsilon _0)}}\frac{1}{v\left(\epsilon \right)}\mathrm{exp}\left[\frac{\left(tT(\epsilon ,\epsilon _0)\right)^2}{4D(\epsilon ,\epsilon _0)}\right].$$
(E6)
At a given time $`t`$ the exponent in the above distribution function is maximal at $`\epsilon =\epsilon _m`$, where $`\epsilon _m`$ is given by the condition $`t=T(\epsilon _m,\epsilon _0)`$, or:

$$t=\underset{\epsilon _m}{\overset{\epsilon _0}{}}\frac{d\epsilon ^{}}{v_0\left(\epsilon ^{}\right)}.$$
(E7)
Expanding the expression (E6) around $`\epsilon =\epsilon _m`$, we have gaussian distribution:

$$P_L(\epsilon ,\epsilon _0,t)\frac{1}{\sigma (\epsilon _m,\epsilon _0)\sqrt{2\pi }}\mathrm{exp}\left[\frac{\left(\epsilon \epsilon _m(\epsilon _0,t)\right)^2}{2\sigma ^2(\epsilon _m,\epsilon _0)}\right],$$
(E8)
with the dispersion:

$$\sigma ^2(\epsilon _m,\epsilon _0)=2v_0^2\left(\epsilon _m\right)D(\epsilon _m,\epsilon _0)=2v^2\left(\epsilon _m\right)\underset{\epsilon _m}{\overset{\epsilon _0}{}}\frac{d\epsilon ^{}}{v_0\left(\epsilon ^{}\right)\mathrm{\Omega }\left(\epsilon ^{}\right)}.$$
(E9)
Note, that according to Eq.(E7), the dispersion is time dependent. If the density of states $`𝒩(\epsilon )`$, and, consequently, $`\mathrm{\Omega }(\epsilon )`$ and $`v_0(\epsilon )`$, are strongly varying functions of $`\epsilon `$, the integrals can be further simplified. In this case we have:

$$\underset{\epsilon _m}{\overset{\epsilon _0}{}}\frac{d\epsilon ^{}}{v\left(\epsilon ^{}\right)}\left[\frac{dv\left(\epsilon _0\right)}{d\epsilon _0}\right]^1\left[\frac{dv\left(\epsilon _m\right)}{d\epsilon _m}\right]^1,$$
(E10)
$$\underset{\epsilon _m}{\overset{\epsilon _0}{}}\frac{d\epsilon ^{}}{v\left(\epsilon ^{}\right)\mathrm{\Omega }\left(\epsilon ^{}\right)}\left[\frac{d\left[v\left(\epsilon _0\right)\mathrm{\Omega }\left(\epsilon _0\right)\right]}{d\epsilon _0}\right]^1\left[\frac{d\left[v\left(\epsilon _m\right)\mathrm{\Omega }\left(\epsilon _m\right)\right]}{d\epsilon _m}\right]^1,$$
(E11)
$$_{\epsilon _m}^{\epsilon _0}\frac{d\epsilon ^{}}{v\left(\epsilon ^{}\right)\mathrm{\Omega }\left(\epsilon ^{}\right)}\left[\frac{dv\left(\epsilon _0\right)\mathrm{\Omega }\left(\epsilon _0\right)}{d\epsilon _0}\right]^1\left[\frac{dv\left(\epsilon _m\right)\mathrm{\Omega }\left(\epsilon _m\right)}{d\epsilon _m}\right]^1.$$
(E12)
The applicability condition for the saddle-point method is:

$$D^3(\epsilon ,\epsilon _0,s_0)\left[\frac{^3}{s_0^3}\underset{\epsilon }{\overset{\epsilon _0}{}}𝑑\epsilon ^{}\frac{\mathrm{\Omega }\left(\epsilon ^{}\right)\mathrm{ln}g(\epsilon ^{},s_0)}{2v\left(\epsilon ^{}\right)}\right]^2.$$
(E13)
Note, that this condition restricts the applicablity of the saddle-point approximation to times, that are not too large. Also, from Eqs.(E5) one can see, that this approximation is invalid when $`\epsilon `$ is close enough to the initial point $`\epsilon _0`$, — here $`\left|s_0\right|`$ becomes large, and it is not possible to expand over $`s_0`$. So, the initial stage of the evolution, when the entire distribution is concentrated near $`\epsilon _0`$, has to be investigated separately too.

## F Calculation of the momenta for large times

Here we present details on the calculation of the momenta of the distribution function for large times and for the exponential density of states. In our calculation we take into account that the density of states decays with decreasing energy.

Let us change the integration variable in (93), $`\rho z=𝒲\left(s\rho e^\rho \right)`$, $`\rho =𝒲\left(s^1ze^z\right)`$. As $`\rho 1`$, we have, introducing $`\rho _s=𝒲\left(s^1\right)`$, $`s=\rho _s^1\mathrm{exp}\left(\rho _s\right)`$:

$`\rho =𝒲\left(\rho _sze^{\rho s+z}\right)\rho _s+z+\mathrm{ln}z,`$

as $`\left|z\right|\left|\rho _s\right|`$ and $`\left|\mathrm{ln}z\right|\left|\rho _s\right|`$. Both inequalities holds within actual region of integration. We can write now:

$$\chi _{nL}(s)=\frac{1}{s}\underset{0}{\overset{\mathrm{}}{}}𝑑z\left[\mathrm{ln}\rho _s\frac{z+\mathrm{ln}z}{\rho _s}\right]^n\frac{}{z}\mathrm{exp}\left[\frac{1}{\overline{\omega }\rho _s^2}z\left(1+\frac{1}{2}z\right)\right].$$
(F1)
Something more about the approximation we have chosen: every extra power of $`z/\rho _s`$ means an extra power of $`\sqrt{\overline{\omega }}`$, which can be neglected. An extra power of $`z`$ means the multiplier $`\sqrt{\overline{\omega }}\rho _s`$. All our considerations are valid if $`\overline{\omega }\rho _s1`$ only. Therefore, we have the reason to omit extra powers of $`z/\rho _s\sqrt{\overline{\omega }}`$ in the following, but to keep some extra powers of $`z^11/\sqrt{\overline{\omega }}\rho _s\sqrt{\overline{\omega }}`$.

Thus we have, for example:

$`\chi _1(s)={\displaystyle \frac{1}{s}}\mathrm{ln}\rho _s+{\displaystyle \frac{1}{s\rho _s}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑z\left(z+\mathrm{ln}z\right){\displaystyle \frac{d}{dz}}\mathrm{exp}\left[{\displaystyle \frac{1}{\overline{\omega }\rho _s^2}}z\left(1+{\displaystyle \frac{1}{2}}z\right)\right]`$ (F2)
$`=`$ $`{\displaystyle \frac{1}{s}}\mathrm{ln}\rho _s+{\displaystyle \frac{\sqrt{\overline{\omega }}}{s}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑y\left(y+{\displaystyle \frac{\mathrm{ln}\left(\sqrt{\overline{\omega }}\rho _s\right)}{\sqrt{\overline{\omega }}\rho _s}}+{\displaystyle \frac{\mathrm{ln}y}{\sqrt{\overline{\omega }}\rho _s}}\right){\displaystyle \frac{d}{dy}}\mathrm{exp}\left({\displaystyle \frac{y^2}{2}}{\displaystyle \frac{y}{\sqrt{\overline{\omega }}\rho _s}}\right),`$ (F3)
where

$`y={\displaystyle \frac{z}{\sqrt{\overline{\omega }}\rho _s}}.`$

Within three leading orders on the parameter $`1/\sqrt{\overline{\omega }}\rho _s`$ we have:

$$\chi _{1L}(s)\frac{\mathrm{ln}\rho _s}{s}\sqrt{\frac{\pi \overline{\omega }}{2}}\frac{1}{s}\frac{\mathrm{ln}\left(\sqrt{\overline{\omega }}\rho _s\right)}{s\rho _s}+\left(1\frac{\mathrm{ln}2\gamma }{2}\right)\frac{1}{s\rho _s}+\sqrt{\frac{\pi }{8\overline{\omega }}}\frac{1}{s\rho _s^2},$$
(F4)
which can be written, using $`\rho _s+\mathrm{ln}\rho _s=\mathrm{ln}\left(1/s\right)`$ also as:

$$\chi _{1L}(s)=\frac{1}{s}\mathrm{ln}\mathrm{ln}\frac{b\sqrt{\overline{\omega }}}{s}\sqrt{\frac{\pi \overline{\omega }}{2}}\frac{1}{s}+\sqrt{\frac{\pi }{8\overline{\omega }}}\frac{1}{s\rho _s^2},b=\sqrt{2}e^{1\gamma /2},$$
(F5)
where $`\gamma `$ is the Euler’s constant. In the time representation we have:

$`\chi _1(t)={\displaystyle _C}{\displaystyle \frac{ds}{2\pi i}}\chi _{1L}(s)e^{st}={\displaystyle _C}{\displaystyle \frac{ds}{2\pi is}}\mathrm{ln}\mathrm{ln}{\displaystyle \frac{b\sqrt{\overline{\omega }}}{s}}e^{st}\sqrt{{\displaystyle \frac{\pi \overline{\omega }}{2}}}+\sqrt{{\displaystyle \frac{\pi }{8\overline{\omega }}}}{\displaystyle _C}{\displaystyle \frac{ds}{2\pi is\mathrm{ln}^2\frac{\sqrt{\overline{\omega }}}{s}}}e^{st}.`$

At large $`t`$, $`b\sqrt{\overline{\omega }}t1`$, the first integral, after the substitution: $`st=y`$, may be expanded as:

$`{\displaystyle _C}{\displaystyle \frac{ds}{2\pi is}}\mathrm{ln}\mathrm{ln}{\displaystyle \frac{b\sqrt{\overline{\omega }}}{s}}e^{st}={\displaystyle _C}{\displaystyle \frac{dy}{2\pi iy}}\mathrm{ln}\left[\mathrm{ln}b\sqrt{\overline{\omega }}t\mathrm{ln}y\right]e^y=\mathrm{ln}\mathrm{ln}b\sqrt{\overline{\omega }}t`$
$`{\displaystyle \frac{1}{\mathrm{ln}b\sqrt{\overline{\omega }}t}}{\displaystyle _C}{\displaystyle \frac{dy}{2\pi iy}}e^y\mathrm{ln}y{\displaystyle \frac{1}{2\mathrm{ln}^2b\sqrt{\overline{\omega }}t}}{\displaystyle _C}{\displaystyle \frac{dy}{2\pi iy}}e^y\mathrm{ln}^2y\mathrm{}`$
$`=`$ $`\mathrm{ln}\mathrm{ln}b\sqrt{\overline{\omega }}t+{\displaystyle \frac{\gamma }{\mathrm{ln}b\sqrt{\overline{\omega }}t}}+{\displaystyle \frac{\pi ^2/6\gamma ^2}{2\mathrm{ln}^2b\sqrt{\overline{\omega }}t}}.`$
As for the second integral, one may restrict oneself with the leading term only. As a result, we have:

$$\chi _1(t)=\mathrm{ln}\mathrm{ln}\stackrel{~}{b}\sqrt{\overline{\omega }}t\sqrt{\frac{\pi \overline{\omega }}{2}}+\sqrt{\frac{\pi }{8\overline{\omega }}}\frac{1}{\mathrm{ln}^2\stackrel{~}{b}\sqrt{\overline{\omega }}t}+\mathrm{},\stackrel{~}{b}=\sqrt{2}e^{1+\gamma /2}.$$
(F6)
Let us consider now the second momentum:

$$\chi _{2L}(s)=\underset{0}{\overset{\mathrm{}}{}}\frac{d\rho }{\rho }\mathrm{ln}^2\rho \varphi _L(\rho ;s)\mathrm{ln}^2\rho .$$
(F7)
Again, we change the integration variable to $`z=𝒲\left(s\rho e^\rho \right)`$, write $`\rho 1`$ as $`\rho =\rho _s+z+\mathrm{ln}z`$, and set the lower limit of the integration to be zero. Restricting ourselves with the lowest order in $`z/\rho _s\sqrt{\overline{\omega }}`$, we set $`\mathrm{ln}\rho =\mathrm{ln}\rho _s+(z+\mathrm{ln}z)/\rho _s`$, and:

$`{\displaystyle \frac{d\rho }{\rho }}\varphi _L={\displaystyle \frac{dz}{s}}{\displaystyle \frac{}{z}}\mathrm{exp}\left[{\displaystyle \frac{1}{\overline{\omega }\rho _s^2}}z\left(1+{\displaystyle \frac{1}{2}}z\right)\right].`$

Taking into account:

$`\chi _{1L}(s)={\displaystyle \frac{1}{s}}\mathrm{ln}\rho _s{\displaystyle \frac{z+\mathrm{ln}z}{\rho _s}},`$

we can write:

$$\chi _{2L}(s)=s\chi _{1L}^2(s)+\left(\frac{z+\mathrm{ln}z}{\rho _s}\right)^2\frac{z+\mathrm{ln}z}{\rho _s}^2.$$
(F8)
Then we have:

$`\left({\displaystyle \frac{z+\mathrm{ln}z}{\rho _s}}\right)^2`$ $`=`$ $`{\displaystyle \frac{1}{s\rho _s^2}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑z(z+\mathrm{ln}z)^2{\displaystyle \frac{d}{dz}}\mathrm{exp}\left[{\displaystyle \frac{1}{\overline{\omega }\rho _s^2}}z\left(1+{\displaystyle \frac{1}{2}}z\right)\right]`$
$`=`$ $`{\displaystyle \frac{\overline{\omega }}{s}}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑y\left(y+{\displaystyle \frac{\mathrm{ln}\left(\sqrt{\overline{\omega }}\rho _s\right)}{\sqrt{\overline{\omega }}\rho _s}}+{\displaystyle \frac{\mathrm{ln}y}{\sqrt{\overline{\omega }}\rho _s}}\right)^2\mathrm{exp}\left({\displaystyle \frac{y^2}{2}}{\displaystyle \frac{y}{\sqrt{\overline{\omega }}\rho _s}}\right).`$
Expanding this expressions in power series on $`1/\sqrt{\overline{\omega }}\rho _s`$, and keeping three leading powers, it may be easily obtained:

$`\left({\displaystyle \frac{z+\mathrm{ln}z}{\rho _s}}\right)^2={\displaystyle \frac{2\overline{\omega }}{s}}+{\displaystyle \frac{\sqrt{2\pi \overline{\omega }}\mathrm{ln}\left(\sqrt{\overline{\omega }}\rho _s\right)}{s\rho _s}}{\displaystyle \frac{\mathrm{ln}2+\gamma }{2}}{\displaystyle \frac{\sqrt{2\pi \overline{\omega }}}{s\rho _s}}+{\displaystyle \frac{\mathrm{ln}^2\left(\sqrt{\overline{\omega }}\rho _s\right)}{s\rho _s^2}}`$
$`2\left(1{\displaystyle \frac{\mathrm{ln}2\gamma }{2}}\right){\displaystyle \frac{\mathrm{ln}\left(\sqrt{\overline{\omega }}\rho _s\right)}{s\rho _s^2}}+\left[{\displaystyle \frac{\pi ^2}{24}}\mathrm{ln}2+\gamma +\left({\displaystyle \frac{\mathrm{ln}2\gamma }{2}}\right)^2\right]{\displaystyle \frac{1}{s\rho _s^2}}.`$
¿From Eq. (F4) we have:

$`{\displaystyle \frac{z+\mathrm{ln}z}{\rho _s}}=\sqrt{{\displaystyle \frac{\pi \overline{\omega }}{2}}}{\displaystyle \frac{1}{s}}+{\displaystyle \frac{\mathrm{ln}\left(\sqrt{\overline{\omega }}\rho _s\right)}{s\rho _s}}\left(1{\displaystyle \frac{\mathrm{ln}2\gamma }{2}}\right){\displaystyle \frac{1}{s\rho _s}}\sqrt{{\displaystyle \frac{\pi }{8\overline{\omega }}}}{\displaystyle \frac{1}{s\rho _s^2}}.`$

So, from Eq. (F8) we obtain:

$$\chi _{2L}(s)=\chi _{1L}^2(s)+\left(2\frac{\pi }{2}\right)+\left(1\mathrm{ln}2\right)\frac{\sqrt{2\pi \overline{\omega }}}{s\rho _s}+\left(\frac{\pi ^2}{24}+\frac{\pi }{2}1\right)\frac{1}{s\rho _s^2}.$$
(F9)
However, we want to calculate central momentum:

$$\mu _2(t)=\chi _2(t)\chi _1^2(t),$$
(F10)
or, in the Laplace representation:

$`\mu _{2L}(s)`$ $`=`$ $`\chi _{2L}(s)\left(\chi _1^2\right)_L(s),`$ (F11)
$`\left(\chi _1^2\right)_L(s)`$ $``$ $`{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑te^{st}\chi _1^2(t)={\displaystyle \underset{i\mathrm{}+0}{\overset{+i\mathrm{}+0}{}}}{\displaystyle \frac{ds_1}{2\pi i}}\chi _{1L}(s_1)\chi _{1L}(ss_1).`$ (F12)
Taking the expression for $`\chi _{1L}`$ from Eq. (F5), we have:

$`\left(\chi _1^2\right)_L(s)={\displaystyle \underset{i\mathrm{}+0}{\overset{+i\mathrm{}+0}{}}}{\displaystyle \frac{ds_1}{2\pi i}}{\displaystyle \frac{\mathrm{ln}\mathrm{ln}\frac{b\sqrt{\overline{\omega }}}{s}\mathrm{ln}\mathrm{ln}\frac{b\sqrt{\overline{\omega }}}{ss_1}}{s_1\left(ss_1\right)}}+{\displaystyle \frac{\sqrt{2\pi \overline{\omega }}}{s}}\mathrm{ln}\mathrm{ln}{\displaystyle \frac{b\sqrt{\overline{\omega }}}{s}}+{\displaystyle \frac{\pi \overline{\omega }}{2s}}`$ (F13)
$`{\displaystyle \underset{i\mathrm{}+0}{\overset{+i\mathrm{}+0}{}}}{\displaystyle \frac{ds_1}{2\pi i}}{\displaystyle \frac{\mathrm{ln}\mathrm{ln}\frac{b\sqrt{\overline{\omega }}}{s_1}}{s_1\left(ss_1\right)\mathrm{ln}^2\frac{b\sqrt{\overline{\omega }}}{ss_1}}}{\displaystyle \frac{\pi }{2s\mathrm{ln}^2\frac{b\sqrt{\overline{\omega }}}{s}}}+\mathrm{}.`$ (F14)
Let us consider the first integral. After changing the integration variable, $`s_1=sx`$, we have:

$`{\displaystyle \frac{1}{s}}\mathrm{ln}^2\mathrm{ln}{\displaystyle \frac{b\sqrt{\overline{\omega }}}{s}}+{\displaystyle \frac{2}{s}}\mathrm{ln}{\displaystyle \frac{b\sqrt{\overline{\omega }}}{s}}{\displaystyle _C}{\displaystyle \frac{dx}{2\pi i}}{\displaystyle \frac{\mathrm{ln}\left(1\mathrm{ln}(1x)/\mathrm{ln}\left(\frac{b\sqrt{\overline{\omega }}}{s}\right)\right)}{x(1x)}}`$
$`+{\displaystyle _c}{\displaystyle \frac{dx}{2\pi i}}{\displaystyle \frac{\mathrm{ln}\left(1\mathrm{ln}x/\mathrm{ln}\left(\frac{b\sqrt{\overline{\omega }}}{s}\right)\right)\mathrm{ln}\left(1\mathrm{ln}(1x)/\mathrm{ln}\left(\frac{b\sqrt{\overline{\omega }}}{s}\right)\right)}{x(1x)}}.`$
Closing the integration contour to the left, one can easily prove, that the first intehral in the above expression is zero. The second one at small $`s`$ may be expanded into powers of $`1/\mathrm{ln}s`$. In the leading order we have:

$`{\displaystyle \frac{1}{s}}\mathrm{ln}^2\mathrm{ln}{\displaystyle \frac{b\sqrt{\overline{\omega }}}{s}}+{\displaystyle \frac{1}{s\mathrm{ln}^2\frac{b\sqrt{\overline{\omega }}}{s}}}{\displaystyle _C}{\displaystyle \frac{dx}{2\pi i}}{\displaystyle \frac{\mathrm{ln}x\mathrm{ln}(1x)}{x(1x)}}.`$

Performing integration by part in the integral, and then deforming the integration contour to make it run along the lower and upper shores of the cut $`(\mathrm{},0)`$ of the integration function:

$`{\displaystyle _C}{\displaystyle \frac{dx}{2\pi i}}{\displaystyle \frac{\mathrm{ln}x\mathrm{ln}(1x)}{x(1x)}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _C}{\displaystyle \frac{dx}{2\pi i}}\mathrm{ln}^2x{\displaystyle \frac{d}{dx}}{\displaystyle \frac{\mathrm{ln}(1x)}{1x}}={\displaystyle \frac{1}{2\pi }}{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑x{\displaystyle \frac{\mathrm{ln}(1+x)1}{(1+x)^2}}Im\mathrm{ln}^2(x+i0)`$
$`=`$ $`{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑x{\displaystyle \frac{\left[\mathrm{ln}(1+x)1\right]\mathrm{ln}x}{(1+x)^2}}={\displaystyle \frac{\pi ^2}{6}}.`$
The evaluation of the second integral in Eq. (F14) gives:

$`2{\displaystyle \frac{\mathrm{ln}\mathrm{ln}\frac{b\sqrt{\overline{\omega }}}{s}}{s\mathrm{ln}^2\frac{b\sqrt{\overline{\omega }}}{s}}}+O\left({\displaystyle \frac{1}{s\mathrm{ln}^4\frac{b\sqrt{\overline{\omega }}}{s}}}\right)`$

As a result, one can write:

$$\left(\chi _1^2\right)_L(s)=s\chi _{1L}^2(s)+\frac{\pi ^2}{6s\rho _s^2}+\mathrm{}$$
(F15)
Thus, one have the second central momentum in the Laplace representation to be:

$$\mu _{2L}(s)=(2\frac{\pi }{2})\frac{\overline{\omega }}{s}+(1\mathrm{ln}2)\frac{\sqrt{2\pi \overline{\omega }}}{s\rho _s}+\left(\frac{\pi }{2}\frac{\pi ^2}{8}1\right)\frac{1}{s\rho _s^2},$$
(F16)
and the time dependence of the second central momentum:

$$\mu _2(t)=\left(2\frac{\pi }{2}\right)\overline{\omega }+\left(1+\mathrm{ln}2\right)\frac{\sqrt{2\pi \overline{\omega }}}{\mathrm{ln}t}.$$
(F17)
Expression for the third central momentum $`\mu _3\left(t\right)`$ (97) may be obtained quite analogously.

## G Asymptotic form of the distribution function at large times

One can write the large times distribution function (92), after the replacement of the integration variable $`sz=𝒲\left(s\rho e^\rho \right)`$as follows:

$$\varphi (\epsilon ,t)=\frac{}{t}\psi (\epsilon ,t)$$
(G1)
$$\psi (\epsilon ,t)=\rho e^\rho \underset{i\mathrm{}+\delta }{\overset{+i\mathrm{}+\delta }{}}\frac{dz}{2\pi iz^2}(1+z)\frac{}{\epsilon }\mathrm{exp}\left[z+t\frac{z}{\rho }e^{z\rho }\frac{1}{\overline{\omega }\rho ^2}z\left(1+\frac{1}{2}z\right)\right].$$
(G2)
One can try to calculate the latter integral within the saddle point approximation. The saddle point equation is:

$`f^{}(z)=1+t{\displaystyle \frac{z+1}{\rho }}e^{z\rho }{\displaystyle \frac{1}{\overline{\omega }\rho ^2}}(z+1)=0.`$

Assuming $`\left|z\right|\overline{\omega }\rho ^2`$ (remind that $`\omega \rho ^21`$), we can neglect the third term in the above equation, and to obtain:

$`z_c+1=𝒲\left({\displaystyle \frac{\rho }{t}}e^{\rho +1}\right).`$

As we shall see later, the actual values of parameters $`\rho `$, $`t`$, corresponding to the body of the distribution function, obeys the inequality $`\rho e^\rho /t1`$, and, therefore, $`z_c1`$. But this leads to the saddle point parameter, $`\left|f^{\prime \prime \prime }(z_c)/\left[f^{\prime \prime }(z_c)\right]^{3/2}\right|=\left|(z_c+3)/(z_c+2)^{3/2}\right|1`$, to be of the order of 1, instead of being much less. Let us note, however, that the function in the exponent in Eq. (G2) may be represented as a sum of two terms, of which the first one, $`z+tz/\rho \mathrm{exp}\left(z\rho \right)`$, ensures the function in the integral to have a maximum at $`z=z_c`$ of the width $`1`$, while the second one, $`\left(\overline{\omega }\rho ^2\right)^1z\left(1+z/2\right)`$, may be supposed to be almost constant within this peak, — it varies essentially at scales $`\mathrm{\Delta }z\sqrt{\overline{\omega }}\rho 1`$ within the actual interval of system’s parameters. Thus, one can replace that second term in the exponent by its value at some $`z_0=z_c+a`$, $`\left|a\right|1`$. After this, the integration may be easily performed, returning back to the initial integration variable $`s`$:

$`\psi (\epsilon ,t)\rho {\displaystyle \frac{}{\rho }}\mathrm{exp}\left\{{\displaystyle \frac{1}{\overline{\omega }\rho ^2}}𝒲\left(s_0\rho e^\rho \right)\left[1+{\displaystyle \frac{1}{2}}𝒲\left(s_0\rho e^\rho \right)\right]\right\}{\displaystyle \underset{i\mathrm{}+\delta }{\overset{+i\mathrm{}+\delta }{}}}{\displaystyle \frac{ds}{2\pi is^2}}e^{st},`$

where the last integral is equal to simply $`t`$, and, taking into account $`z_01`$, $`s_0=z_0/\rho \mathrm{exp}\left(z_0\rho \right)c/t`$, where $`c1`$ will be obtained later. Performing differentiation with respect to $`t`$, see Eq. (G2), and neglecting $`z_0/\overline{\omega }\rho ^21`$, we obtain finally:

$$\varphi (\epsilon ,t)\rho \frac{}{\rho }\mathrm{exp}\left\{\frac{1}{\overline{\omega }\rho ^2}𝒲\left(\frac{c}{t}\rho e^\rho \right)\left[1+\frac{1}{2}𝒲\left(\frac{c}{t}\rho e^\rho \right)\right]\right\}.$$
(G3)
To obtain $`c`$, let us calculate $`\chi _1(t)=\overline{\epsilon }(t)`$ with the distribution function (G3). We have:

$`\chi _1(t)={\displaystyle \frac{d\rho }{\rho }\varphi (\rho ,t)\mathrm{ln}\rho }={\displaystyle 𝑑z\mathrm{ln}\rho (z)\frac{}{z}\mathrm{exp}\left[\frac{1}{\overline{\omega }\rho ^2(z)}z\left(1+\frac{1}{2}z\right)\right]}`$

$`\rho (z)=𝒲\left({\displaystyle \frac{t}{c}}ze^z\right),`$

where the replacement of the integration variable $`\rho z=𝒲\left(c\rho e^\rho /t\right)`$ was done. Assuming $`\rho 1`$, one can write, using $`𝒲\left(x\right)x\mathrm{ln}x`$ as $`\left|x\right|1`$, that $`\rho (z)𝒲\left(t/c\right)+z+\mathrm{ln}z\rho _t+z+\mathrm{ln}z`$ (the condition $`z\sqrt{\overline{\omega }}\rho _t\rho _t=𝒲\left(t/c\right)`$ was used). In the leading order on $`\sqrt{\overline{\omega }}`$, one can replace $`\rho (z)`$ within the exponent with $`\rho _t`$, and $`\mathrm{ln}\rho (z)\mathrm{ln}\rho _t+(z+\mathrm{ln}z)/\rho _t`$. The limits of integration by $`z`$ are to be set $`(0,+\mathrm{})`$. Then we have:

$$\chi _1(t)=\mathrm{ln}\rho _t\underset{0}{\overset{\mathrm{}}{}}𝑑z\frac{z+\mathrm{ln}z}{\rho _t}\frac{}{z}\mathrm{exp}\left[\frac{1}{\overline{\omega }\rho _t^2}z\left(1+\frac{1}{2}z\right)\right].$$
(G4)
This latter integral, analogously to one in Eq. (F), may be evaluated as an expansion on small parameter $`1/\sqrt{\overline{\omega }}\rho `$:

$$\chi _1\left(t\right)=\mathrm{ln}\mathrm{ln}\left(\sqrt{2}e^{1\gamma /2}\frac{\sqrt{\overline{\omega }}t}{c}\right)\sqrt{\frac{\pi \overline{\omega }}{2}}+\mathrm{}$$
(G5)
Comparing Eqs.(G5) and (95), one can see, that $`c=\mathrm{exp}\left(\gamma \right)`$.