Generalized phase-space kinetic and diffusion equations for classical and dispersive transport

We formulate and solve a physically-based, phase space kinetic equation for transport in the presence of trapping. Trapping is incorporated through a waiting time distribution function. From the phase-space analysis, we obtain a generalized diffusion equation in configuration space. We analyse the impact of the waiting time distribution, and give examples that lead to dispersive or non-dispersive transport. With an appropriate choice of the waiting time distribution, our model is related to fractional diffusion in the sense that fractional equations can be obtained in the limit of long times. Finally, we demonstrate the application of this theory to disordered semiconductors.


Introduction
The link between theory and experiments measuring transport properties in either gases or condensed matter is provided by an advection-diffusion equation or Fokker-Planck equation for number density n t r ( , ) [1]. These configuration-space theories are only valid in the hydrodynamic regime of smooth spatial gradients [2]. It is therefore fortunate that physical systems such as crystalline condensed matter and gaseous media rapidly approach this regime. In these systems, collisions are effectively instantaneous, and memory of the initial condition is lost after only a few collisions. Any large gradients are quickly smoothed. For such 'classical' transport, the corresponding diffusion equation is a conventional partial differential equation of first order in t and second order in r. Such systems are characterized by a diffusive regime wherein the mean square displacement grows linearly, i.e. 〈 〉 − ∼ x x t 2 2 . In contrast, there is an increasing body of literature demonstrating 'anomalous' diffusion, wherein the mean square displacement is nonlinear, i.e. 〈 〉 − ∼ γ x x t 2 2 where γ ≠ 1 [3]. To mention just a few, subdiffusion (γ < 1) occurs in the transport of charge carriers in disordered semiconductors [4,5], and the movement of lipids and proteins of cell membranes [6]; whereas Lévy flight (γ > 1) has been observed in the diffusion of ultracold atoms in an optical lattice [7,8].
Subdiffusion is fundamentally slower than regular Brownian diffusion. It arises from repeated trapping and de-trapping of the transported species, in which periods of 'classical' transport are interrupted by (potentially long) periods of immobilization. Consequently, the memory of the initial condition may persist for long times, and large gradients are not necessarily smoothed out. Such effects are often accounted for by replacing the conventional time derivative with a fractional derivative [9][10][11][12][13], which in turn accounts for the memory effects. However, the spatial gradient terms are left intact, and thus the weak gradient assumption is therefore still implicit. This is an apparent contradiction that challenges the validity of the fractional diffusion equations. Until now, this issue has only been addressed in an ad hoc matter, through solution of an exactly solvable model kinetic equation in phase space [14]. The phase space system does not require the assumption of weak gradients, so this questionable assumption is avoided. Trapping effects were incorporated in the 'collision term' in a way which, although consistent with other contemporary kinetic equations, is nevertheless ad hoc in nature, and therefore warrants further scrutiny [14].
This paper resolves the ad hoc nature of the previous approach by revisiting the phase space formulation. We develop from first principles a new physically-based, exactly solvable model kinetic equation, whose solution leads, in the weak gradient limit, to a generalized diffusion equation. This equation provides a general description of transport in the presence of trapping. In the appropriate limit of instantaneous de-trapping, it reduces to the classical diffusion equation. Conversely, in the case of traps that possess a divergent waiting time, dispersive (subdiffusive) transport is obtained, similarly to the fractional diffusion equations commonly studied [15][16][17][18].

Procedure
Kinetic theory aims at finding the charge carrier phase space distribution function f t r v ( , , ), from which the quantities of physical interest follow as integrals over velocity space, e.g., the number density 3 . In the classical kinetic theory of Boltzmann, collisions are assumed to be instantaneous in time [19], but here collisions times may be nonzero. We quantify the collision times according to a probability density function, which we call the 'waiting time distribution' ϕ τ ( ). A generalized kinetic equation for the phase-space distribution function f t r v ( , , ) is derived, from which we obtain an exact expression for the number density n t r ( , ). In the weak gradient limit, n is shown to satisfy a generalized diffusion equation, valid for any ϕ τ ( ), which can be applied to either classical or trap-limited transport. 3. The phase space model Figure 1 outlines the phase space trapping and de-trapping picture schematically. Suppose that at time t charge carriers moving freely are scattered out of a phase space element centred on r v ( , ) at a constant rate ν into localized or trapped states, i.e., the rate of loss is -νf t r v ( , , ). They remain trapped for a range of various possible times τ, as determined by the waiting time distribution, defined such that ϕ τ τ ( ) d is the probability of de-trapping between times τ and τ τ + d after the trapping event. Thus the rate at which particles are re-entering the phase-space is the distribution function for the de-trapped particles, assumed to be a Maxwellian at medium temperature T. Here, m is the particle mass (or effective mass), and k B is the Boltzmann constant. The corresponding free particle function is thus where a is the external force per unit mass. Note that only the second (de-trapping) term on the right hand side is convolved with ϕ, corresponding to a delayed release from the localized states. In contrast, scattering into traps takes place without any delay, and the first term on the right hand side is not convolved. We note that the classical BGK equation [20], well known in gas and crystalline semiconductor transport studies, is regained when collisions are instantaneous. Instantaneous collisions are recovered by setting the waiting time distribution to be the Dirac delta, i.e. ϕ τ δ τ = ( ) ( ). We now define the operator t t and write equation (2) in the equivalent form Integration of (4) over all velocities yields the equation of continuity, is the free particle flux. A further integration over r yields is the total free particle number. Since this implies (2) in infinite space, with the initial condition,

The solution of equation
corresponding to release of N (0) particles from the origin of coordinates with a Maxwellian distribution of velocities with an initial temperature T 0 , can be obtained exactly through Fourier and Laplace transformation in space and time respectively, following the same mathematical procedure as [14] and [21]. The transformed number density then follows is the Laplace transform of the waiting time distribution. Inversion of the Laplace transform of equation (8) could be carried out, if desired, through the usual contour integral, which would be evaluated using the residue theorem in terms of the singularities of n p k ( , ) , which are given by the zeroes of the denominator of (8), i.e., However, here the interest is confined to the asymptotic, weak gradient, long time regime, which is determined by the small | | k solutions of (10). Since in this case 2 4 of the plasma dispersion function [22] may be used. Proceeding in this way, the solution of (10), valid to second order in | | k , is found to be where a and k are column vectors, I is the unit matrix, and : denotes a double contraction over tensor indices 1 . A similar result follows from Laplace-Fourier transformation of the generalized diffusion equation t d with ∂ t specified by equation (3). The singularities of n p k ( , ) in this case are found from d T which, to be consistent with equation (12), requires the drift velocity and diffusion tensor to be given by It is clear that the effect of trapping enters equation (14) only through the operator ∂ t , while the spatial gradient terms are determined by free carrier transport. Moreover, the free carrier transport coefficients are unaltered by trapping, e.g., equation (16) provides exactly the same expressions that one obtains from the classical (non-trapping) BGK model kinetic equation [20,21]. At this point we note that we can obtain the same result in a more direct way, by simply assuming Fickʼs law for the free carriers, d and substituting into the right hand side of equation (5), furnishing equation (14) once more. From this perspective, the preceding phase space analysis may be taken as confirming the validity of Fickʼs law even in the presence of trapping and de-trapping.
To summarise this section, we have developed a microscopic phase-space kinetic theory including trapping. This theory was physically motivated and did not require the ad hoc introduction of 'fractional' terms to incorporate memory effects, as was done previously [14]. By solution of this microscopic theory, we have identified the regime of validity for the macroscopic Fickʼs law, and consequently, justify the generalized diffusion equation (14). In the following section, we proceed to apply this model to a time-of-flight experiment.

Solution of the generalized diffusion equation
Consider now a slab of material of thickness L between two plane-parallel electrodes, the normal direction defining the z-axis of a system of coordinates. Assuming all spatial dependence to be in this direction only, equations (14) and (3) together yield the generalized diffusion equation in one dimension For an impulse initial condition, δ = − n z N z z ( , 0) (0) ( ) 0 , and perfectly absorbing boundaries, = = n t n L t (0, ) ( , ) 0, we solve this by firstly taking the Laplace transform and then applying the Poisson summation theorem as outlined in [23], to obtain  To this point the discussion is quite general, but to go further, we must specify ϕ t ( ).

Waiting time distributions ϕ(t ) with a finite first moment-'classical' transport
To study the impact of the waiting time distribution on the solution, we initially considered several simple ϕ t ( ) functions. Time of flight transients were calculated using equation (21), numerically inverting the Laplace transform. Dimensionless results are shown in figure 2, where time is scaled to the transit time = t L v / tr d through a sample of thickness L, and concentration is scaled to the initial concentration N (0). In this system of units, the normalized drift velocity is equal to 1. We set the normalized diffusion coefficient to The simplest choice for the waiting time distribution is the Dirac delta, ϕ δ = t t ( ) ( ). In this case, de-trapping occurs instantaneously and consequently has no impact. The system reduces to the standard advection-diffusion equation.
Another choice is an exponential distribution, ϕ τ = , where again τ is the first moment of the waiting time distribution function. As can be seen from figure 2, the first moment of ϕ controls the long-time behaviour. Higher moments can only influence the behaviour at shorter times, as is demonstrated by the differences between the cases with an exponential distribution and those with the Laplace domain series expansion.

Relating the waiting time distribution to a density of trapped states
Rather than assuming an ad hoc waiting time distribution, it may be useful to calculate it from a more fundamental physical model. In what follows, we give a specific example for how this might be achieved. We consider a semiconductor with traps that form a density of localized states below the band gap. The release times ϕ t ( ) are determined by the distribution of these traps in energy space.
To describe this semiconductor, we use a multiple trapping model with a uniform capture (trapping) cross-section [24] for charge carriers. We define the density of localized states to be g(E), where < E 0 is the energy relative to the conduction band. If the rate of escape from a trap at energy E is proportional to where ν 0 is a frequency characterizing the rate of escape from traps. The density of states g(E) can be measured experimentally [25,26]. In this case we assume an exponential distribution, which occurs in organic and inorganic materials [27,28].
where T c is a characteristic temperature that describes the width of the density of states. Equation (22)  where γ · · ( , ) is the lower incomplete Gamma function 2 , and α = T T / c . This distribution appears in the literature of the multiple trapping model, for example, equation (9) of [29]. This distribution is normalized, and it has a divergent first moment, sufficient to describe dispersive transport [15,30,31]. Figure 3 shows typical time of flight transients based on equation (19) together with the waiting time distribution equation (23). Initially, all carriers are assumed to be untrapped. At short times the profiles are classical, with a transition to dispersive behaviour at longer times as the carriers begin to enter the trap states. In the dispersive regime, the sum of slopes is −2, exactly as expected for an exponential density of states [4]. We note that a different form of dispersive transport may arise from a different density of states, and that our model can be readily adapted by evaluating equation (22) with the appropriate g(E) function.
The parameters used for figure 3 were chosen as follows. The attempt to escape frequency ν 0 is extremely fast (e.g. ∼10 12 Hz [24]), so we consider only the situation where ν ≫ t 1 tr 0 .