Third-order transport coefficients for localised and delocalised charged-particle transport

We derive third-order transport coefficients of skewness for a phase-space kinetic model that considers the processes of scattering collisions, trapping, detrapping and recombination losses. The resulting expression for the skewness tensor provides an extension to Fick’s law which is in turn applied to yield a corresponding generalised advection-diffusion-skewness equation. A physical interpretation of trap-induced skewness is presented and used to describe an observed negative skewness due to traps. A relationship between skewness, diffusion, mobility and temperature is formed by analogy with Einstein’s relation. Fractional transport is explored and its effects on the flux transport coefficients are also outlined.

Very little data regarding third-order transport coefficients (the skewness tensor) can be found in the literature. This is understandable, since they have not been included in the interpretations of traditional swarm experiments. There is, however, a growing interest regarding these transport coefficients, partially due to estimations that third-order transport coefficients could be measured in the present or near future 1,2 . It is also considered that third-order transport coefficients would be very useful, in combination with transport coefficients of a lower order, for determination of cross section sets, by means of inverse swarm procedure 1,2 . Third-order transport coefficients are also required for the conversion of the hydrodynamic transport coefficients into transport data measured in steady state Townsend and arrival time spectra experiments 3,4 . The skewness tensor can also be employed in fluid models of discharges, by pairing a generalised diffusion equation, which includes the contributions of the third-order transport coefficients, with Poisson's equation. This could be particularly important for discharges where ions play an important role 5 , or in situations where the hydrodynamic approximation is at the limit of applicability, as in the presence of sources and sinks of particles or in the close vicinity of physical boundaries.
In this study, we are concerned with the form of the skewness tensor for charged-particle transport in the presence of trapped (localised) states. In particular, we are interested in the scenario where transport is dispersive. Dispersive transport is characterised by a mean squared displacement that increases sublinearly with time 6 . Due to this non-integer power-law dependence, we refer to dispersive transport as fractional transport throughout this study. Some examples of fractional transport include the trapping of charge carriers in local imperfections in semiconductors [7][8][9][10][11] and both electron [12][13][14] and positronium [15][16][17] trapping in bubble states within liquids. Third-order transport coefficients are expected to be more sensitive to the influence of non-conservative collisions than those of lower order, suggesting that the presence of such trapped states would significantly influence the skewness tensor. Indeed, skewness and other higher order transport coefficients are used to characterise fractional transport in a variety of contexts, including transport in biological cells [18][19][20][21] . Consider also Fig. 1, which plots the solution of the Caputo fractional advection-diffusion equation, a common model for fractional transport 6 . This solution exhibits a large skewness in comparison to the accompanying Gaussian solution of the corresponding classical advection-diffusion equation.
In the following, we describe charged particle transport using a full phase-space kinetic model as defined by a generalised Boltzmann equation with a corresponding trapping and detrapping operator. In our previous papers [22][23][24] , we introduced and studied such a generalised Boltzmann equation, deriving lower-order transport coefficients (up to diffusion) and generalisations of the Einstein relation. We will extend the results of these papers to determine the skewness tensor. Calculations of the skewness tensor for the Boltzmann equation have been performed previously by a number of authors 2,25-28 . We will use these earlier studies to confirm the structure of the skewness tensor and to benchmark our results in the trap-free case.
In Sec. 2 of this study, we outline a general phase-space kinetic model 23 for charged-particle transport via localised and delocalised states. This model is capable of describing both normal and fractional transport. We follow in Sec. 3 with a derivation of the flux transport coefficients for this model up to third order. Sec. 4 explores the structure of these transport coefficients and their symmetries under parity transformation. The transport coefficients are used to extend Fick's law, which leads to a generalised advection-diffusion-skewness equation, presented in Sec. 5. In this section, we also provide a physical interpretation of trap-induced skewness. By analogy with Einstein's relation, Sec. 6 provides a relation between skewness, diffusion, mobility and temperature. Sec. 7 looks at the case of fractional transport and its effects on the flux transport coefficients. Finally, Sec. 8 lists conclusions along with possible avenues for future work.

Phase-space kinetic model
We previously reported 22-24 the development of a phase-space kinetic model wherein charged particles scatter due to collisions, enter and leave traps and undergo recombination. In this model, free particles are described by a phase-space distribution function f(t, r, v), defined by the generalised Boltzmann equation where E is the applied electric field and particles have charge e, mass m and number density n(t, r) ≡ ∫ dvf (t, r, v).
Here, collisions, trapping and free particle recombination occur at the constant frequencies ν coll , ν trap and ν loss (free) , respectively. For collisions, the Bhatnagar-Gross-Krook (BGK) collision operator 29 has been used, which scatters particles isotropically according to a Maxwellian velocity distribution of background temperature T coll . We define the Maxwellian velocity distribution of temperature T as where k B is the Boltzmann constant. Similarly, we use the BGK-type operator introduced by Philippa et al. 22 to describe the processes of trapping and detrapping. This operator specifies that particles leave traps with a Maxwellian distribution of velocities of temperature T detrap after a delay that is governed by the distribution of trapping times φ(t). That is, trapping events are described mathematically as delayed scattering events. This distribution appears in Eq. (1) through the effective waiting time distribution

Transport coefficients to third order
By integrating the generalised Boltzmann equation (1) throughout all of velocity space, we find the equation of continuity for the number density n(t, r): trap loss (free) where the particle flux is In the weak-gradient hydrodynamic regime, physical quantities can be written as an infinite series of spatial gradients of the number density n(t, r) 30,31 . In the case of the flux Γ(t, r), such a density gradient expansion provides a generalisation of Fick's law: 2 where W is the drift velocity vector, D is the rank-2 diffusion tensor and Q is the rank-3 skewness tensor. To determine these flux transport coefficients it is simply a matter of writing the solution of the generalised Boltzmann equation (1) itself as a density gradient expansion and then evaluating the flux using Eq. (5), resulting in the transport coefficients Substituting the density gradient expansion of f (t, r, v) into the Boltzmann equation (1) and equating coefficients of spatial gradients, as done in Sec. IV of ref. 23 , gives the following coefficients coll trap (1) where a Fourier transform has been performed in velocity space, . As in ref. 23 , we have used the density gradient expansion of the concentration of particles leaving traps: (1) ( 2) 2 the coefficients of which are related to the flux transport coefficients through Applying this time average to unity results in an implicit definition for the initial coefficient R: Thus, for every trapping time distribution φ(t) there corresponds a value of R. Some values are tabulated in Appendix A of ref. 23 .
Proceeding to evaluate Eqs (8-10) for the transport coefficients, we find where e 1 , e 2 and e 3 are standard orthonormal basis vectors and we have introduced the effective frequency and temperature: We confirm that when there are no traps present, ν trap = 0, the transport coefficients agree with those of the BGK collision model, previously found by Robson 26 :

Structure and Symmetry of Transport Coefficients
If we align the basis vector e 3 parallel to the applied electric field E, the transport coefficients (19-21) take on the known tensor structure 2,25,28,30,31 :  where a, b ∈ {x, y, z}. Here, the drift velocity is defined by the speed eff the diffusion coefficient is defined by two components perpendicular and parallel to the field trap trap 2 eff and the skewness is defined by the three independent components Although this is the case in general, there are situations where the skewness can be defined using fewer than three components. Indeed, this is the case for the BGK model as studied by Robson 26 where the skewness given by Eq. (26) is defined using only the components Q 1 and Q 3 , with Q 2 = 0. The component Q 2 vanishes in this case due to the simple Maxwellian source term used to describe scattered particles. For Q 2 to arise, it is necessary that this source term has some spatial dependence, as occurs for our model through the concentration of particles leaving traps, ϕ * t t r ( ) n( , ), and its density gradient expansion (14). Lastly, we also confirm that the symmetry of transport coefficients with respect to the parity transformation E → −E depends on the parity of the order of each transport coefficient 25,32 :

Generalised Advection-diffusion-skewness Equation
Using the density gradient expansion (6) for the flux Γ(t, r) up to second spatial order in conjunction with the continuity equation (4) results in the generalised advection-diffusion-skewness equation valid in the weak-gradient hydrodynamic regime. In Cartesian coordinates (x, y, z) with the electric field E aligned in the z-direction, the transport coefficients take the form of Eqs (27)(28)(29)(30)(31) where the skewness manifests as components perpendicular and parallel to the applied field 2,5,28 : zzz which in terms of the independent components (35-37) are Written in full, the perpendicular and parallel skewnesses are  where terms present due to trapping have been grouped separately and the lower-order transport coefficients (32)(33)(34) have been used to simplify. An alternative form of the skewness tensor that makes use of these components explicitly is zab where a, b ∈ {x, y, z}. This form was used by Robson 26 when expressing the BGK model skewness (26) and is valid only when the skewness is triple-contracted with a symmetric tensor, as occurs in the advection-diffusion-skewness equation (41).
To provide some physical intuition regarding the perpendicular and parallel skewness coefficients, Q ⊥ and Q , we solve the advection-diffusion-skewness equation (42) for an impulse initial condition and perform contour plots of the resulting pulse in Fig. 2. Figure 2(a) considers the case of no skewness, = = ⊥ Q Q 0, and displays the expected Gaussian solution with elliptical contours due to anisotropic diffusion. Figure 2(b) and (c) consider positive perpendicular and parallel skewnesses, respectively. In both cases, it can be seen that skewness introduces asymmetry in the pulse in the direction of the field. In general, positive skewness can be seen to reduce the spread of particles behind the pulse, while enhancing the spread toward the front of the pulse. In Fig. 2(b) for positive perpendicular skewness, this change in particle spread primarily occurs transverse to the field, resulting in a vaguely triangular pulse profile. In Fig. 2(c) for positive parallel skewness, this change in particle spread occurs longitudinally which, in the language of statistics, results in a distribution with positive skew.
In our previous paper 23 , we interpreted the trap-induced anisotropic diffusion present in Eq. (34) as a consequence of the physical separation between trapped particles and free particles moving with the field. In a similar fashion, we can interpret the trap-induced skewness present in the perpendicular and parallel skewness coefficients (47) and (48). To achieve this, we plot the skewness against the detrapping temperature T detrap for various mean trapping times in Fig. 3. The resulting plots are linear with gradients that characterise of the type of skewness caused by traps. That is, positive or negative gradients correspond respectively to positive or negative trap-based skewness.
When the mean trapping time is zero, the gradients in Fig. 3 are positive and traps cause positive skewness. This is to be expected as, in this case, trapping and detrapping simply act as an elastic scattering process with a positive skewness akin to Eq. (26) for the BGK collision model. As the mean trapping time increases, the nature of the skewness caused by traps changes, ultimately becoming negative for the parameters considered in Fig. 3. As illustrated in Fig. 2, negative skewness corresponds to an increased spread of particles behind the pulse. We interpret the increased spread here as being due to particles returning from traps. This interpretation implies that the skewness coefficients could become overall negative if particles remain trapped for a sufficient length of time before returning with a sufficiently large temperature. Indeed, these are the conditions for which the skewness coefficients become negative in Fig. 3.
This phenomenon of negative skewness has been observed previously by Petrović et al. 5 in the calculation of the perpendicular skewness of electrons in methane. Only collisions were considered in this study and so trapping is evidently not a necessary condition for negative skewness to occur. However, it should be emphasised that the skewness is strictly positive when collisions are described by the simple BGK collision operator, as is seen in Eq. (26).

Relating skewness, mobility and temperature
The classical Einstein relation between diffusion, mobility and temperature is 33   . To perform these plots, we choose a trapping frequency of ν trap /ν coll = 1/9, while (b) also requires that we specify a drift velocity W, which we choose such that mW 2 /k B T coll = 1/4. The gradients in (b) are of smaller magnitude than (a) due to the greater dependence of the parallel skewness (48) on the drift speed W as compared to the perpendicular skewness (47). Thus, as the drift speed decreases, the plots in (b) coincide with those in (a). When detrapping is instantaneous, τ = 0, the skewness gradients are positive, implying that the skewness caused by traps is also positive. As the mean trapping time τ increases, the skewness gradients decrease, becoming negative and implying a corresponding negative skewness due to traps. The limiting case of an infinite mean trapping time, τ → ∞, corresponds to fractional transport, which is the subject of Sec. 7. We observe from this figure that the skewness coefficients become overall negative when particles leave traps with a sufficiently large temperature T detrap after a sufficiently long amount of time τ. This observation coincides with the illustration of skewness in Fig. 2 where negative skewness is characterised by an increased particle spread behind the pulse, which we attribute here to particles returning from traps.  Koutselos 34 has presented a similar relationship between the skewness tensor and lower-order transport coefficients for the case of the classical Boltzmann equation.

The case of Fractional Transport
For the phase-space kinetic model described by Eq. (1), fractional transport can occur when the distribution of trapping times has a heavy power-law tail of the form 23 (1 ) Note that, as transport here is dispersive in nature, the mean trapping time diverges:  Note that transport coefficients are now independent of the specific choice of the trapping time distribution φ(t), so long as the condition (55) for fractional transport is satisfied.
Knowledge of the skewness coefficient (59), as well as other higher-order transport coefficients, is useful for characterising fractional transport. For example, Norregaard et al. 18 use higher-order moments to analyse the motion of biological particles.

Conclusion
We have explored the transport coefficients of a phase-space kinetic model (1) for both localised and delocalised transport. In particular, we have considered up to the third-order transport coefficient of skewness bfQ, which takes the form of a rank-3 tensor. The structure of the skewness tensor and its symmetry under parity transformation was found to be in agreement with previous studies. These transport coefficients provide an extension to Fick's law, Eq. (6), which we used to form a generalised advection-diffusion-skewness equation (41) with a non-local time operator. We observed trap-induced negative skewness and provided a corresponding physical interpretation. In addition, by analogy with Einstein's relation, the skewness was related to the mobility and temperature through Eq. (54). Lastly, the form of the transport coefficients for the particular case of fractional transport were outlined in Eqs (57-59).
There exist a number of possibilities for future work. The focus of this paper was on constant transport coefficients that define the flux in the hydrodynamic regime as the density gradient expansion (6). Transient transport coefficients and transport coefficients of the bulk were not considered. Ref. 23 outlines an analytical solution of the kinetic model (1) that could be used to compute such transport coefficients through time-varying velocity and spatial moments of the phase-space distribution function f (t, r, v). Another extension to this work could be to explore what consequences energy-dependent collision, trapping and recombination frequencies have on the skewness. Such a generalisation for Eq. (1) was the focus of ref. 24 . This would allow for the derivation of a skewness analogue of Einstein's relation that would also take into account the field dependence of mobility 24 . This may also shed light on the recent results of Petrović et al. 5 , that suggest a correlation between the energy-dependent phenomenon of negative differential conductivity and skewness.
Lastly, it is important to note that the extension to Fick's law described in this paper is only useful when an electric field is present. Without an applied field, the drift velocity, skewness and all other odd-ordered transport coefficients would vanish. If we wish to extend Fick's law in such a situation, we must also consider the kurtosis coefficient, the next even-ordered transport coefficient beyond diffusion. The kurtosis can be found in a straightforward fashion from the rank-3 tensorial coefficient f (3) (v) in the density gradient expansion (7) of the phase-space distribution function f (t, r, v), in the same way drift velocity, diffusion and skewness were found using Eqs (8)(9)(10). Data availability statement. No datasets were generated or analysed during the current study.