Particle propagation and electron transport in gases

In this review, we detail the commonality of mathematical intuitions that underlie three numerical methods used for the quantitative description of electron swarms propagating in a gas under the effect of externally applied electric and/or magnetic fields. These methods can be linked to the integral transport equation, following a common thread much better known in the theory of neutron transport than in the theory of electron transport. First, we discuss the exact solution of the electron transport problem using Monte Carlo (MC) simulations. In reality we will progress much further, showing the interpretative role that the diagrams used in quantum theory and quantum field theory can play in the development of MC. Then, we present two methods, the Monte Carlo Flux and the Propagator method, which have been developed at this moment. The first one is based on a modified MC method, while the second shows the advantage of explicitly applying the mathematical idea of propagator to the transport problem.


Introduction
The goal of kinetic theory is the modelling of a gas (or plasma, or any system made up of a large number of particles) by a distribution function in the particle phase space.This phase space includes macroscopic variables, i.e. the position in physical space, but also microscopic variables, which describe the "state" of the particles.In this respect, several textbooks and review papers have been devoted to describe the mathematical foundations of kinetic theory of gases and plasmas [1][2][3][4][5][6].The need to create mathematical models of ionized gases at low temperatures has also been felt for many decades.In fact, the applications of these systems are numerous and qualitatively important and range from the detection of ionizing radiation to the production of new materials, to medicine, and the emission of light [7].It is not possible to build a reliable model of a low-temperature ionized gas without a detailed description of electron transport and kinetics.For this reason, a gigantic literature was born and developed, including thousands of works, on the development of numerical calculation methods that allow us to calculate the quantities that describe these phenomena [8][9][10].
Since the density of electrons in a weakly ionized gas is much lower than that of the neutral background gas, it is usually assumed that they do not collide altogether or with other charged particles.Their behavior is governed by externally applied fields (electric or magnetic) and by the various collisions that they undergo with the molecules or atoms of the neutral gas.Under these assumptions, it follows that the motion of electrons in the phase space can be described by a one-particle distribution function f (r, v, t) under the linear electron Boltzmann equation (EBE) [11,12]: where r, v, and t are the configuration space coordinates, velocity space coordinates, and time, respectively, e is the elementary charge, m is the electron mass, E is the electric field, and J[f ] is the collision term that takes into account electron collisions with the neutral background gas atoms or molecules.For the category of problems in question, this equation which is a non-linear equation becomes linear because in the product of distribution functions in the right-hand-side of the equation, one of the two functions describes the neutral gas particles.Even with this simplification, the equation presents considerable difficulties and this is why the literature relating to its solution is so extensive [8,10,13,14].
For many years, and still effectively when appropriate, an approximation known as the two-term approximation has been employed to solve this equation [11,[15][16][17][18].It consists in simplifying the electron velocity distribution function (EVDF) by developing it into an expansion in spherical harmonics truncated at first order.This approximate methodology is quick and simple to implement numerically and has been available for several years also into open access programs, which allow one to calculate the isotropic and the first anisotropic components of the velocity distribution function [17,19,20].
Over the years, it has emerged that the two-term approximation is not sufficient in a series of problems especially related to plasmas confined inside reactors, whose geometry must be considered [21][22][23][24][25][26].For this reason and with the help of the ever-improving performance of computers, more direct calculation methods have been developed.These allow the distribution function of Equation ( 1) to be determined if necessary, taking into account all the independent variables involved.
The first of these methods, still one of the most used today, is the Monte Carlo (MC) method which directly simulates the trajectory of each electron in a large ensemble using random numbers [27][28][29][30].This method is a variation of the one originally developed in the 1940s for calculating neutron trajectories in the simulation of nuclear fission systems [31].Nowadays, several codes are available as open source for MC simulations of electrons in gases [32][33][34][35].Despite the accuracy of MC simulations for calculation of chemical rate coefficients and electron transport parameters and the increase in computational resources, this method is still computationally expensive (especially when coupled with self-consistent description of weakly-ionized gases) [36].
Recently, other deterministic methods for numerical solutions of the EBE have emerged [37][38][39].For a long time, these methods were essentially unaffordable because of the high computational cost of the integral in the Boltzmann collision operator, but they are now becoming more and more competitive.The methods that will be exposed in this review share the role played in their formulation and practice by mathematical entities called Green functions or propagators [40].These, basically, are operators that allow to calculate the evolution of a system, as long as it is described by linear equations, by breaking it down into the contributions that come from every small part of the system into every other small part of the system.Green functions, on which there is a gigantic literature, place the methods considered here in the broader perspective of mathematical physics.Furthermore, awareness of the common conceptual basis between these methods produces advantages both for those who develop them, for those who employ them and for those who teach them.In the following sections we will develop this point of view.
To summarize, the numerical methods developed for solving the EBE are an important part of mathematical physics and theoretical physics, and any research in this sector deserves to be considered in the broader context.On the occasion of the 150 th anniversary of the original formulation of the EBE, a review article on methods for calculating the kinetics of electrons based on the EBE was published by Boyle and co-authors [10].The present review has been formulated as a complement to that work.In fact, we focus here on methods for numerical solutions of the EBE that do not employ an expansion of the velocity distribution function in spherical harmonics.
Moreover, the present review updates the survey written by Segur and co-authors [13] with a description of recent efforts in development of numerical methods used to solve the EBE.
The review is structured as follows.Section 2 describes the principles of MC simulations of electrons and the connection of MC methods with the electron transport equation.Section 3 presents the mathematical principles behind the Monte Carlo Flux method, that is a hybrid stochastic-deterministic algorithm to solve the EBE.The advantages and disadvantages of this method are also described.Finally, in Section 4, the Propagator method is described and its application for electron swarm and plasmas are highlighted, as well as the numerical method.To conclude, we show that awareness of the common conceptual basis between these methods produces advantages both for those who develop them, for those who employ them, and for those who teach them.In the following sections, we will develop this point of view.

Principle
The MC method is still today the most used for the description of electron swarms especially in situations of complex geometry or when, in the context of a simulation of a plasma device, electrons are interacting with space charge or with the walls of the device [9,28,29,33,36,[41][42][43][44][45].
The method is based on the description of the movement of a large number of mathematical objects representing electrons under the effect of forces due to electric and magnetic fields, while collisions with neutral plasma particles are introduced by means of random times between each collisions and the next one.
It is very simple to create a MC simulation [27]; it is sufficient to organize a vector of objects representing the electrons or, in a more traditional setting, numerical vectors for the individual dynamic quantities of each electron, i.e. x, y, z, v x , v y , v z .Once a small, appropriately chosen time step has been introduced, it is sufficient to take into account the effect of the forces acting on each electron due to its position and speed.At this point, it is sufficient to apply the kinetic theory of gases to determine the probability that a given collision process extracted from the set of cross sections can occur.The effect of this collision process on any single electron is produced explicitly, for example by modifying the velocity components to take into account an elastic or inelastic process.
Alternatively, each electron is associated with the time to its next collision, this is updated during the calculation, and when this is reached the speed of the electron is changed and a new time to the next collision is calculated.the version of the method most often used for precision calculations is the one called null collisions MC.The method was first introduced by Skullerud [46] for the use in plasma, although a similar one, the artificial isotope method [31], had been previously used for the treatment of complex geometries in nuclear reactors.The basic idea is simple and powerful and consists in introducing a collision process, fictitious, in such a way that its frequency added to that of the other processes produces a collision frequency that does not depend on space, time or speed.This way it is extremely simple to calculate the time to the next process, only once after any collision event, using the equation [47]: where t c is the time to next collision, ν is the above mentioned constant total momentum transfer collision frequency (including the null collision frequency), and η is a random number from a uniform distribution in (0, 1].To compensate for the fictitious process that has been added, this last is simply considered equal to the others, but when it occurs, the speed of the electron does not change.This method also allows the effect of the non-zero temperature of the gaseous medium on the electrons to be introduced: it is clear in fact that the collision frequency between a neutral particle and an electron depends on the speed of the particle.This effect is particularly important if ions are propagated instead of electrons [48,49].A reasonable maximum of the collision frequency is then taken, and at the time of the collision a fraction of the events are eliminated to account for the Maxwell distribution of gas particle velocities.

Monte Carlo method as formal solution of the electron transport problem
The MC method, with the prescription summarized above, is able to solve the problem of the transport of a swarm in an exact way, given that no numerical parameter is introduced that needs to be optimized.For this reason it has sometimes been considered a kind of numerical experiment based on an analogical approach, on the direct use of the laws of physics, something different from the solution of transport equations.This idea is in contrast with the very principles of the initial development of the MC method, at the time of which it was perfectly clear to the developers that the method can be derived by formal procedures starting from the integral form of the transport equation [31]; there is therefore no opposition between an analogical approach and an approach based on the solution of the integral equation, especially since the two methods ultimately lead to the development of the same code.
In the case of applying the method to swarms of charged particles, there are some formalities to introduce in the mathematical proofs, given that charged particles do not propagate in straight lines unlike photons and neutrons.However, these technical aspects have been addressed in several publications [2,50].In this regard, an illuminating way of showing how the continuous operators of the EBE become events in the MC method is to use the well-known integral equation [50] satisfied by the time evolution operator U written as the exponential of the sum of two operators At and Jt where t is time: very well known in quantum field theory, but mathematically true in general; Now if exp(−At) is the operator which propagates a particle for the time t under the action of inertia and electric force, while −J is the integral operator representing the collision events, The equation above shows that the last collision event randomly breaks the propagation into two parts, or stages, one of which, U (t ′ ), can still include other collision events in a hierarchical structure.
This structure is resumed in a natural way using a diagramming technique that was mentioned in the work [27] developed in more detail in [51].
To show this, we introduce, for the first time explicitly, the Green function or propagator between two positions in phase space after a time shift ∆t = t ′ − t that is represented here as G(r ′ , r, v ′ , v, t ′ , t) which actually is a generalized function, since it may include the Dirac's delta δ(r The intuitive aspect of traditional MC methods, derives from the fact that in the absence of collisions the Green function of two arguments final v ′ , r ′ , t and initial v, r is equal to 0 unless there is an electron trajectory, determined by inertia and by the electric and possibly magnetic fields, which smoothly joins the two arguments; in which case the Green function has the aforementioned form of the Dirac generalized function.
Calculating this propagator without collisions is therefore equivalent to solving the equations of motion for a material point with a charge-to-mass ratio corresponding to an electron, under the effect of fields.
When collisions of electrons with atoms and molecules are added to the description, the point probability packets are dispersed, which means that the initial Dirac distribution becomes a Dirac distribution multiplied by a value less than 1, added to a traditional function.At this stage, it is no longer possible to use a description based only on deterministic trajectories.Now, we can introduce here a notation well known in quantum field theory [52] where G, the Green function of the equation, is the exact propagator, while we use g to represent the Green function where the positive part of the collision operator J is removed, so that particles propagate according to the convective operator while they "decay" and disappear according to the exponential factor exp(−νt).The result is the expansion: The expression above is quite complex, and its subsequent orders are increasingly cumbersome, but it can be conveniently summarized by diagrams as shown in Figure 1.In the figure, once again using a notation borrowed from quantum field theory [53], Figure 1: Representation of a MC calculation as a perturbative expansion of the exact propagator using a diagrammatic technique, explained in the text.A MC trajectory with a certain number of collisions n c is a contribution to the calculation of the diagram with n c dots.After [51].
an arrow represents a Green function, a dot represents the effect of integration and of the positive part of the collision operator J, a thick arrow represents the exact Green function G (i.e. the Green function of the full transport equation).g, the thin arrow, represents not only collisionless particle propagation, but includes exponential decay, this is essential for the full G, the thick arrow, to conserve the normalization of f .If the null-collision method is used, the operator represented by the dot will contain a term not affecting f and a term with cJ, c < 1.
The factor 1/2 which precedes the second term in Equation ( 6) takes into account the fact that a collision at t 1 followed by a collision at t 2 is equivalent to the same collisions in opposite order.The term of n-th order will therefore be preceded by a factor 1/(n!).If the perturbation expansion is written using the null collision method, then a factor of the type exp(−νt)[νt] n /(n!) will appear in front of each term in the expansion [49,51]; the details of the derivation are sketched in the next section.These factors constitute a Poisson distribution, which can be generated using the formula for collision times (2).In this way, the formula normally introduced on an empirical basis finds a mathematical motivation.
These diagrams show that the null-collision MC method is in fact directly related to a specific form of the transport equation [47].Indeed, it is possible to formally state that, given a set of MC simulations of a given case study, the subset that includes exactly n-collisions corresponds to a stochastic calculation of the term of n-th order in a perturbation expansion of the solution of the equation transport.
Sometimes, one can be misled by the fact that the MC method introduces statistical fluctuations of the quantities then felt in the transport equation, such as in particular the kinetic distribution of velocities.This seems to show that it has a higher element of describing reality than these equations.In reality the statistical fluctuations of the MC method are not physical and simply have to be reduced as much as possible by accumulation of events, leading to statistical convergence towards the exact solution.
The method always arrives at this for an infinite number of events, when the solved equation is linear.We take advantage of this aspect to note something that is also essential in the other sections of this review, namely the fact that swarm problems are generally described by linear equations, given that the particles of the swarm are diluted to the point of not interacting.This is a basic requirement to apply the concept of the Green function.
These considerations allow us to place the MC method, normally considered a technical tool in the field of plasma simulation, in a broader context which is that of methods for theoretical physics.A very appropriate position given that it is formally capable of solving the Boltzmann transport equation, and not only as a numerical experiment, from an analogical point of view, as the terminology often used in publications seems to suggest [2,6].
This physical-mathematical aspect of the method can be disregarded in its most basic implementations and it can be developed precisely having in mind a direct image of the physics of the transport process.

Principle
In the previous section, we noted that the traditional MC method corresponds in any case to a formal calculation, i.e. the calculation of the Green functions, or electron propagators.These functions are generalized functions (i.e.Dirac delta functions), with a diffuse component which becomes the majority as time increases, precisely due to collisions [51].The computed Green functions in each case depend on continuously varying arguments, which represent initial and final position and velocity components.
Since the beginning of the '90s, modifications of the method have been experimented, in which instead of calculating propagators at the level of description of the electrons for the entire simulated time, the formalism of Markov processes has been employed to separate the diffusion time scale from the collision scale [54].Among these, the Monte Carlo Flux method (MCF) stands out in particular, which is a hybrid method between the deterministic method of the propagators and the stochastic MC method [37].The fundamental idea of this method is the separation between diffusion and collisional time scales.In fact, due to the small mass of the electron compared with that of atoms and molecules, in some very common regimes, hundreds of collisions are needed to observe a significant average dispersion of the initial kinetic energy.It is therefore convenient to describe the diffusive phase by means of averaged quantities obtained previously from calculations describing the effect of the collisions on the single electrons.To make this calculation practical, it is advisable to carry out a coarse-grained preliminary operation by introducing a lattice of suitable geometry into the configuration space.Once this is done, diffusive kinetics over long times can be described using the Master Equation or Kolmogorov equation [54,55] (i.e. the balance of probability over time).This corresponds to the Markov process in continuous time but in a discretized space.
The schematics of MCF is shown in Figure 2, where the three main steps of the numerical methods are highlighted (in order from left to right).The transition frequencies p ij (∆t) which appear as coefficients depending on two indices and on time in the equation are precisely the quantities which can be calculated by means of a collection of traditional MC simulations [37].In each of these collections of simulations, electrons are uniformly distributed within a cell of the lattice in the configuration space and a simulation of the duration ∆t corresponding to a few collisions is therefore performed.
At the end of the simulation, the electrons found in each of the cells estimate the transition frequency between the initial cell and the final cell.Running this collection of preliminary simulations can take significant time, as it is necessary to break down for each individual cell.However, the Markov matrix that is calculated in this way can be used as a Green function in coarse-grained space and for a fixed time difference, so that it is possible to perform different simulations, for different times as long as they are multiples of ∆t, and for several initial distributions, all without the need to run the preliminary MC simulations again.Alternatively, the stationary distribution can be calculated by looking for the stationary vector with respect to the p ij matrix.All of these calculations are completely deterministic and pretty fast.

Mathematical derivation of Monte Carlo Flux
In this subsection, the MCF method is derived starting from a continuous generalization and proceeding with its discretized form, obtained under the multigroup approximation [56].This subsection is based on Chapter 3 of the PhD thesis of one of the authors [57].
The system under consideration is assumed to be homogeneous, such that spatial terms in the EBE can be neglected.Under these assumptions, the EBE is where f (v, t) is the EVDF, a(t) is the electron acceleration due to an external electric field, and J[f ] is the collision term, that in its stochastic form is Moreover, the following usual normalization condition applies for the transition probabilities with ν(v) the total electron impact collision frequency.
A possible technique for the solution of Equation ( 7) is based on the path-integral formulation, as described by Rees [58].In [58], the rigorous path-integral theory has been initially developed for transport studies in semiconductors.The same theory has been extended by Kumar [59] and Longo [47] for the motion of charged particles in a gas in the presence of an external field.The use of this formulation has a twofold advantage, that is; (i) the null-collision technique [46] is directly derived from the transport problem and (ii) the solution of Equation ( 7) is obtained by iterative applications of a continuous generalization of a Markov operator on the initial EVDF.In particular, from Equation ( 8) and Equation ( 9), it is possible to verify that Equation ( 7) can be rewritten as where and ν max is the maximum collision frequency defined such that ν max ≥ ν(v) for all values of the electron speed v.Note that the first term on the right-hand-side of Equation (11) represents the gain of electrons having velocity v at time t due to collisions with the background gas, and the term (ν max − ν(v)) represents "null-collision" events per unit time that leave the EVDF unchanged.The definition of ν max is also widely used in MC simulations, since it is at the basis of the null-collision technique introduced by Skullerud [46].For stationary conditions the solution of Equation ( 10) is The method for evaluating the right-hand-side of Equation ( 12) is based on the following iterative technique [58].First, an intermediate function Second, the distribution The advantage of this method is that the n-th iterate f (n) (v) is equivalent to the distribution obtained after a time interval n/ν max .Furthermore, the stationary EVDF can be obtained in the limit As described by Nanbu [60], the same iterative approach can be applied when the temporal evolution of the EVDF is sought after.In particular, for sufficiently high values of ν max , let ∆t = (ν max ) −1 be the time step for EVDF evolution.Then, at first order in ∆t, the EVDF can be expressed as or equivalently If ∆t is sufficiently low, an iterative solution of Equation ( 10) describes also the temporal evolution of the EVDF, starting from an initial distribution f (0) , at time t = 0.
To summarize, it is useful to represent the aforementioned equations in a schematic form.In particular, the integration on the right-hand-side of Equation ( 13) can be represented as the action of an operator Ŝ on f (n−1) (v) and that on the right-hand-side of Equation ( 14) by the action of an operator P on M (n−1) (v).A complete iteration is therefore described by the action of the combined operator M = Ŝ P on f (n−1) (v).
Moreover, the time evolution of an initial function f (0) (v) after a time interval n/ν max is determined by { M } n .Different numerical techniques can be employed for calculations of the aforementioned operator.In the MCF method, a grid is defined for calculations of transition probabilities for the electron motion between velocity space cells.In this discretized form, M is also termed as the transport matrix (or Markov matrix).
The following subsections are devoted to the description of this discretization and the numerical procedure that is applied for calculations of M .

Discretization of the Electron Velocity Distribution Function
We describe here the numerical method that is used to implement the MCF.A more detailed comparison of MCF against other methods, such as conventional MC, twoterm Boltzmann solver, and multi-term Boltzmann solver, is described in [37,61,62].
In this subsection, a uniform external electric field E is considered, such that Hence the electron velocity v = (v x , v y , v z ) can be represented in polar coordinates (ϵ, θ, ϕ), where ϵ is the electron kinetic energy (i.e.ϵ = (1/2)mv 2 ), θ is the angle between v and the z-direction, and ϕ is the azimuthal angle around the z-axis.Under the assumption of isotropic scattering around the z-axis, the velocity space can be described by only two variables (ϵ, θ), instead of three.Moreover, under this assumption, the EVDF is uniform in ϕ (i.e.f (v) = f (ϵ, cos θ)).Hence the velocity space (ϵ, cos θ) is partitioned into cells C i,j , with 1 ≤ i ≤ I the index for the energy component and 1 ≤ j ≤ J the index for the angular component.The calculation range is assumed as 0 ≤ ϵ ≤ ϵ max and −1 ≤ cos θ ≤ 1, where ϵ max is chosen such that we can neglect electron diffusive fluxes to energies higher than ϵ max .In this way, the energy and angular bin sizes are defined as ∆ϵ = ϵ max /I and |∆(cos θ)| = 2/J.
Under these assumptions, the EVDF can be represented in its discretized form as a column vector f of size IJ × 1 as where each element of the vector f has size I × 1 and represents the total number of particles having energies between 0 and ϵ max and cos θ = cos θ j , with j = 1, . . ., J.In particular, the j-th component (n j ) is written as where n i,j is the total number of electrons in the C i,j cell and it is computed as where N elec is the total number of particles in the MC simulation and In other words, after discretization of the energy and angular components, each simulated particle contributes to the EVDF as a Kronecker delta function [28].From the EVDF, the l-th order Legendre polynomial coefficients f l (ϵ i ) can be calculated as where A l (ϵ i ) = (2l + 1)/(n e ∆ϵ √ ϵ i ) is a normalization factor that depends on the total electron population n e , such that f ⊤ (ϵ i ) is the transpose vector for the EVDF computed at ϵ = ϵ i , thus having size 1 × J, and L l is a column vector of size J × 1 written as P l (cos θ 2 ) . . .
with P l (cos θ j ) the l-th order Legendre polynomial calculated at cos θ j .Note that the Electron Energy Distribution Function (EEDF) comes from Equation ( 22) as the zero-th order (isotropic) component (f 0 ) of the EVDF expansion.Hence f 0 can be calculated with the following simplified expression: meaning that the EEDF can be retrieved directly from binning in energy of the total number of simulated particles.In the following, the basic idea underlying MCF simulations is covered, that is the calculation of the temporal evolution of the EVDF, given an initial distribution at time t = 0.

Calculation of the transport matrix
In the MCF method, the temporal evolution of the EVDF is obtained by a matrix-vector operation, given an initial distribution function at the time t = 0.The EVDF at time t + ∆t is computed as [37] f where M (∆t) is the transport matrix of size IJ × IJ.Note that the matrix M is a discrete version of the continuous operator M defined in Section 3.2.
More generally, if the reduced electric field and the gas composition are fixed for any t > 0, the transport matrix remains unchanged and the EVDF can be found with an iterative procedure as with n the iteration number, and n∆t the time-step associated with the n-th iteration.
An important consequence of equation ( 27) is that the evolution of the system after ∆t is determined only by the state of the system at a time t and it is not affected by the previous history.This is known as Markov property and allows one to rewrite the linear EBE (Equation ( 7)) as a simple Markov chain consisting of the system of linear equations in (27).This property is typically satisfied if t coll ≪ ∆t < t SS , that is ∆t should be much longer than the time interval between two successive collisions t coll , but also shorter than the time t SS for the distribution function to reach steady-state.In particular, collisions are essential for the randomization of the particle velocities and trajectories.Through this randomization, electron history is erased and the evolution of the system depends only on the current state, not on past states.Hence, as shown in [61,62], the choice of an appropriate ∆t is fundamental to ensure the applicability of Equation ( 27).
The form of the transport matrix is where m j,j ′ are sub-matrices of size I ×I.The total number of sub-matrices in Equation ( 28) is determined by the angular discretization.For example, if J cells are defined for the angular component, then all the J × J possible combinations are included in M .
Each sub-matrix is written as: where p i,i ′ |j,j ′ (∆t) is the transition probability for electrons moving from cell C i,j to C i ′ ,j ′ within the time interval ∆t.In the MCF method, short MC simulations are performed for calculations of transition probabilities between all velocity space cells.
Hence, transition probabilities in Equation ( 29) are calculated as where n i ′ ,j ′ (∆t) is the total number of electrons moving from cell C i,j to C i ′ ,j ′ within the time interval ∆t and n i,j (0) is the total number of electrons in the initial cell C i,j at time t = 0. Note that the time interval ∆t depends on physical parameters such as the electric field E and the gas number density N .Empirically, it becomes shorter at higher E and N , as the electron collision frequency increases.The choice of an optimal ∆t is important for an optimization of the CPU time in MC simulations, as well as for accurate calculations of transition probabilities [61].

Advantages and disadvantages
The main advantages of the MCF method are the following [57,61,62]: • Since the MCF is used to calculate the transition frequency and not directly the kinetic distribution, the MCF results have uniform statistical fluctuations for all the regions of the distribution, in particular also the tail, which instead may be statistically inaccessible to the traditional method given the low number of electrons described.
• The matrix-based approach does not require a series expansion of the EVDF and it is easier to implement than a multi-term Boltzmann solver.In perspective, the MCF method can be combined with efficient algorithms for matrix operations and GPU acceleration.
• As a difference with respect to other variance reduction techniques, mainly based on variable mathematical weights for the simulated particles, the MCF method addresses the fundamental problem of the very large ratio, amounting to several orders of magnitude, between the relaxation time of the distribution and the intercollision time.Hence, the stochastic part of MC simulations is limited to a small time interval, that is typically orders of magnitude lower than the steady-state time for the electron energy distribution function (EEDF).
However, in spite of the advantages of this method in terms of both computational cost and accuracy with respect to a conventional MC method, a number of disadvantages of MCF should be considered when choosing an appropriate computational method for plasma applications.For example, • The increasing size of the transport matrix lengthens the computational time and this is also a limit for practical use of the MCF method.This is a problem, for example, for an extension to the configuration space, since additional calculations of transition probabilities of electrons moving between different cells in the spatial coordinates are needed.Nevertheless, matrices in this case are largely sparse.
Hence, the use of efficient algorithms for sparse matrix calculations could help to reduce the large memory that is required.
• The method presented here can only deal with calculations of flux transport parameters.However, bulk parameters are needed when comparing results of calculated electron transport coefficient with swarm measurements, especially at high E/N [67].An extension of the MCF method to the configuration space will provide calculations of the aforementioned parameters as well.
• The transport matrix includes the effect of both field and collisional events.This has practical limitations for calculations of transition probabilities in the presence of a time-varying electric field evolving in timescales comparable with the energy relaxation time.
Some of the limitations of the MCF method can be overcome by the Propagator method that is described in the following Section, where the the field acceleration is separated from the collision part by defining different propagator matrices similar to the transport matrices used in MCF.This makes the Propagator method also applicable to AC electric fields.However, the the Propagator method is limited to low ∆t such that the electron out flow from a cell does not exceed the total number of electrons in that cell [68].Furthermore, it is important to highlight that, as a difference with respect to the Propagator method and the Convection scheme [69], the MCF uses MC simulations for calculations of transition probabilities between velocity space cells.The use of MC simulations is advantageous for its simpler implementation with respect to efficient convective schemes, but it is computationally expensive for large domain simulations.

Principle
Propagator method (PM) is a numerical scheme to calculate the EVDF (or EEDF) on the basis of the EBE.The PM shares many common concepts with the MCF method.The PM calculation is performed with cells defined by partitioning of phase space (v, r), velocity space and real space, into small sections in which electrons There are various cell configurations and experimental models for the PM calculations depending on the target properties of electron swarms.Let us overview efforts of calculations which have ever been performed using the PM.

Classification of cell configurations and expressions of electron motion
We may classify the types of the PM configurations on the basis of the partitioning of velocity space for cells and the expressions of the electron motion in quantification of the intercellular electron transition.
For the former point, the cells can be defined, for example, for every ∆v x , ∆v y , and ∆v z in Cartesian coordinate system (v x , v y , v z ) (Figures 3a and 3b), or for every ∆v, ∆θ, and ∆ϕ in polar coordinate system (v, θ, ϕ) (Figures 3c and 3d), where v z = v cos θ, v x = v sin θ cos ϕ, and v y = v sin θ sin ϕ.As well, division for every ∆ϵ instead of ∆v is also possible.They are referred to as Cartesian-v, polar-v, and polar-ϵ configurations of the cells for velocity space in the following sections.

Collision propagator
The collision propagator represents the change of electron velocity due to collision and scattering.The collisions are categorized into mainly four types.In the simplest case, the following processes are assumed for the electrons undergoing collision at an energy ϵ ′ .

Elastic collision
The electrons are scattered isotropically without loss of energy.They are redistributed to the destination cells having the same energy ϵ = ϵ ′ in proportion to the solid angle of the destination cell subtended at the origin of velocity space (Figure 4a).

Excitation collision
The electrons lose excitation energy ϵ exc and are redistributed to the lower-energy cells of ϵ = ϵ ′ − ϵ exc (Figure 4b).

Ionization collision
The electrons lose ionization energy ϵ ion and the residual energy is shared by the primary and secondary electrons.The electrons, which are doubled, are redistributed to the lower-energy cells of ϵ ≤ ϵ ′ −ϵ ion with relevantly given ratios under the law of energy conservation (Figure 4c).

Electron attachment
The electrons captured by gas molecules disappear from velocity space (Figure 4d).
The number of electrons undergoing collisions of the kth kind in a cell during ∆t, n e,coll,k , is evaluated as n e,coll,k = [N q k (v)v∆t]n e,cell , where N is the gas molecule number density, and q k (v) is the electron collision cross section of the kth kind of collisional process as a function of the electron speed v.The portions n e,coll,k are subtracted from n e,cell of the source cell, and distributed to the destination cells corresponding to the electron velocity after scattering.∆t must be short enough so that the probability of multiple collisions during ∆t can be neglected, and to satisfy k n e,coll,k < n e,cell .

Models of velocity-space under uniform electric fields
The PM calculation for the temporal development of the EVDF of an electron swarm under a uniform electric field E = (0, 0, E z ) is the simplest self-consistent model.It assumes an electron swarm in boundary-free real space and the position of each electron is ignored.The E field may be dc, ac (typically rf: radio frequency), or even impulse fields.
The EBE for this model is written as The spatial electron distribution has already been integrated and its profile is not cared here.
The EVDF in this model can be assumed to be in a rotational symmetry for the azimuthal direction ϕ around the v z -axis.Thus, the EVDF can be represented as a two-variable distribution as f (v ∥ , v ⊥ , t) or f (v, θ, t).Here, v ∥ = v z = v cos θ and Some of the earliest PM efforts to obtain f (v, t) adopted the Cartesian-v configuration with ∆v ∥ and ∆v ⊥ under dc [70,71] and rf [72] E fields.Such a cell configuration simplifies the calculation of the electron acceleration because it is a parallel shift to the +v ∥ direction in velocity space.Each cell has a boundary in contact with its downstream neighbour, and the boundary is normal to a.The probability of electron transition from a cell can be calculated easily in both Eulerian and Lagrangian manners.
The Cartesian-v configuration has also been adopted in calculation of f (v, z, t) extended to a parallel-plane electrodes model including one-dimensional real space z [73][74][75].In the one-dimensional real-space model, there are different treatments for simultaneous changes of v and z in an electron flight, which is explained in the next subsection.
On the other hand, polar-v or polar-ϵ configuration with ∆v or ∆ϵ, and ∆θ or ∆(cos θ) makes the redistribution of electrons after collisions simple when isotropic electron scattering is assumed.The ratio of the scattered electrons that a destination cell receives is proportional to the solid angle of the cell subtended at the origin v = 0 of velocity space in case of isotropic scattering, and the sold angle Ω of a cell is given from dΩ = 2πd(cos θ) = 2π sin θdθ.Such a polar-ϵ configuration has been adopted not only for f (v, t) in velocity space in pulsed Townsend (PT) mode [68,76] but also for the steady-state Townsend (SST) mode [77][78][79], for the time-of-flight (TOF) mode [68,80,81], and rf [82] and impulse [83] fields as mentioned afterward.

Models of real space between parallel-plane electrodes
Phase space to consider the EVDF in one-dimensional real space between parallel-plane electrodes is required to be three-variable for example as (v ∥ , v ⊥ , z) or (v, θ, z).In the point of required memory capacity for practical calculations, it is beneficial that the range of z is finite between the electrodes.The computational array for the cells becomes three-dimensional, and cells are prepared for every ∆v ∥ , ∆v ⊥ , and ∆z [73][74][75], or for every ∆ϵ, θ, and ∆z [77].The same collision propagator as used in the velocity-space model can be applied to this mode because the scattering changes only v of electrons and their positions are unchanged at that moment.The acceleration propagator becomes to involve the spatial displacement because v and z change at the same time in electron free flight under an electric field.
In some early PM efforts using cells defined for every ∆v ∥ , ∆v ⊥ , and ∆z [73][74][75], the destination cells were chosen under a concept of Lagrangian cell.The choice of destination cells and evaluation of the electron redistribution ratios for them are easy in the Cartesian-v configuration because the geometrical shape of cells is simple.This treatment also arrows flexible change of E not only for dc E but also for rf E [74,75] or position-dependent E [73,74].
On the other hand, treatment of electron displacement in an Eulerian manner was also attempted [77] adopting restrictions for the cell configuration and the intercellular electron motion.The cells are defined for every ∆ϵ, ∆θ, and ∆z to satisfy ∆ϵ = eE∆z.
The number of electrons flowing out of a source cell is evaluated by integrating the electron flux passing through the downstream cell boundary in velocity space.For the electrons which move from the source cell defined at an ϵ region to destination cells at the ϵ ± ∆ϵ regions, the spatial displacement ±∆z accompanies, respectively.The changes of ϵ and z are strictly bound under the law of conservation of energy.A demonstration was made under a SST condition by superposing a temporal development of an isolated electron swarm starting from the cathode with an initial energy ϵ = 0 until almost all electrons are absorbed by the anode.Under the restrictions to guarantee the energy conservation, position-dependent profiles of spatial relaxation processes of mean electron energy and drift velocity agreed with results of a MC simulation.Especially, the rising position of ionization coefficient was successfully reproduced at z = ϵ ion /(eE), where ϵ ion is the ionization threshold, by suppressing the numerical diffusion which may occur among the destination cells in the Lagrangian manner.
The boundary condition at the electrodes is also a point of discussion.The simplest model assumes perfect absorption, where the electrons reaching the electrodes disappear.
Elastic or inelastic reflection at the electrodes can be considered by relevantly choosing v after reaching either of the electrodes.The effect of secondary electron emission [74] can be considered as a practical condition as well by supplying new electrons to the cells near the electrode on which primary electrons impinge.

Models of boundary-free real space in steady-state Townsend condition
When SST condition is assumed for one-dimensional electron flow in z direction under uniform electric field E = (0, 0, −E), partitioning of the z position can be omitted [78] in the PM calculation for the EVDF by utilizing the assumptions that the normalized EVDF is identical irrespective of z and that the electron number density varies exponentially with z as n e (z) = n e (0) exp(αz), where α is the effective ionization coefficient.The EVDF is available from the PM calculation performed with only cells for velocity space.Such a technique is to reduce the required memory capacity.
If the EVDF in a slab z 0 ≤ z ≤ z 0 + ∆z is in equilibrium, the electron flow at z = z 0 + ∆z is exp(α∆z) times as large as that at z.The PM calculation considers only f (v, θ) in the slab.A polar-ϵ configuration was adopted and ∆z was set to be ∆ϵ/(eE) to take account of the conservation of energy [78].Electrons in the slab may flow out of the slab, but other electrons flow into the slab from the opposite side to keep the electron population in the slab unchanged.The outflow can be evaluated directly from the EVDF in the slab in the same way as done between parallel-plane electrodes.The inflow is evaluated using the assumption of exponential spatial growth.Let n f,out and n b,out be the numbers of electrons flowing out of the slab forward (v z > 0) at z = z 0 + ∆t and backward (v z < 0) at z = z 0 , respectively, and n f,in and n b,in be those flowing into the slab forward at z = z 0 and backward at z = z 0 + ∆z, respectively, during a time step ∆t.Their relations are The value of α is unknown at this moment.However, with n ion and n att respectively being the electron increase and decrease due to ionization and attachment during ∆t, the value of exp(α∆z) is obtained from a quadratic equation derived from the electron conservation in the slab.They satisfy This treatment is equivalent to the assumption ∂/∂z = α and ∂/∂t = 0 in two-term approximation of the EBE analysis for the SST condition [84].α is obtained from The positive sign is taken in usual SST condition so that α = 0 when n ion − n att = 0.By a relaxation of the EVDF starting from initial EVDF and α, the EVDF in equilibrium under the SST condition is obtained.
Another solution of the quadratic formula corresponding to the negative sign is understood to represent properties of electrons in backward diffusion.In this case, α no longer represents the exponential growth of n e (z) by ionization, but it represents the decay of n e (z) toward the −z direction in the upstream region from the electron source [79,85].A PM calculation with such α [79] demonstrated that the EVDF obtained in the upstream region represents the properties of the missing electrons [86] forming the decay of n e (z) in front of the absorbing anode [73,77].

Models of boundary-free real-space time-of-flight condition
The spatial electron distribution p(z, t) = v f (z, v, t)dv can be composed by superposing the kth-order Hermite functions H k (Z) exp(−Z 2 ) up to a sufficient order n [80]: where H k (Z) is the kth-order Hermite polynomial derived sequentially from where δ ij is the Kronecker delta.The Hermite functions are suitable to expand or compose Gaussian-like distributions having boundary condition lim z→±∞ p(z, t) = 0, and it is thought that p(z, t) of electron swarms eventually tends to a Gaussian.Values of w k (t) are determined by the spatial moments of the electron swarm up to the kth order as shown later.
Let m k (v, t) be the kth-order spatial moment distribution function with respect to the z direction, and it is defined as A series of moment equations in a hierarchy are derived from the EBE by integrating each term with weight z k over z [68,80,81], and m k (v, t) can be calculated without the partitioning of z: Here, m 0 (v, t) is identical to f (v, t) representing the electron number density at v. m 1 (v, t) gives the center of mass G(v, t) of the electrons having v as G(v, t) = m 1 (v, t)/m 0 (v, t), and the center of mass G(t) of the whole electron swarm is given as . From these quantities we can obtain higher-order (nth-order, n ≥ 3) diffusion coefficients D Ln as well [81]; e.g., The same technique to compose spatial electron distribution has been applied not only to the longitudinal direction [81] but also to the transverse direction [89].
The higher-order transverse diffusion coefficients D Tn are also available from the simultaneous moment equations up to the nth order with respect to the direction perpendicular to the E field [89].

Models of velocity space and real space under uniform electric and magnetic fields
The PM has been applied to calculation of the EVDF in equilibrium under dc crossed electric and magnetic fields, E × B fields, assuming E ⊥ B as E = (0, 0, −E) (E > 0) and B = (0, B, 0) (B > 0) [90].The EVDF is no longer axi-symmetric under the E × B fields, thus it is required to have three variables to represent the EVDF.A polar-ϵ configuration modified for three-variable velocity space (v, θ, ϕ) was chosen in the practical calculation for convenience in the treatment of isotropic scattering after collisions [90].The memory array for the EVDF becomes three-dimensional as the cells were prepared for every ∆ϵ, ∆θ, and ∆ϕ, where v = v 1 ϵ/ϵ 1eV , v 1 is the electron speed associated with ϵ 1eV = 1 eV, v x = v sin θ cos ϕ, v y = v cos θ, and v z = v sin θ sin ϕ.A cell may have at most six boundaries facing to the ±ϵ, ±θ, and ±ϕ directions.The EBE for the EVDF in velocity space becomes The treatment for the collision operator is unchanged even for the EVDF in the E × B fields.On the other hand, the acceleration −(e/m)(E+v×B) is velocity-dependent and electron motion in velocity space is rotational around an axis (v x , v z ) = (v E×B , v E ) = (E/B, 0).This rotation axis is parallel to the v z axis but off-centered (off-origin).This arose a complexity in the preparation of the acceleration propagator.The probability of the intercellular electron transition is calculated in the Eulerian manner common with that applied to two-variable velocity space, i.e. by integration of the outflow flux at the downstream cell boundary.However, the cell boundaries facing to the 'downstream' depends on the acceleration direction determined by v. Further more, the acceleration direction may change even in a cell boundary.The acceleration propagator to deal with the rotational electron flow was prepared choosing the downstream cell boundaries carefully, and the order of application of the acceleration propagator to the cells are also arranged for fast convergence of the EVDF.
In an example of computation of the EVDF under E ×B fields with 1000×45×720 cells for (ϵ, θ, ϕ) in double precision (8 bytes per cell) [90], the required memory capacity was roughly several GiB.In addition to the array to store the numbers of electrons in the cells, those for associated cell properties such as cell volume and areas of the six cell boundaries are also needed.However, the memory capacity required for this model is still in an available range for recent workstations.
The EBE under the E × B fields can be extended also to the spatial moment equations under the E × B fields with respect to the x, and y, and z directions [91] in the same manner as Equation ( 41) performed in the boundary-free one-dimensional TOF model.The center-of-mass drift velocity vector W r = (W x , W y , W z ) and the direction-dependent diffusion coefficients D x , D y , and D z were calculated by the PM for the simultaneous spatial moment equations up to the second order.Seven sets of cells were used to store the zeroth-, first-, and second-order moments with respect to the x, y, and z directions, and m z2 (v, t), where m 0 (v, t) is common for the three directions.The collision and acceleration propagators are unchanged.Temporal development of the moments can be calculated by applying the propagators iteratively.On the other hand, in case only the equilibrium values of W x , W y , W z , D x , D y , and D z are needed, the relaxation of each moment can be achieved stepwise from the zeroth order to the second independently for x, y, and z, because higher-order moments depend on only the lower-order ones defined for the same direction via the drift term.This feature allows us to save the computational load for the relaxation process and the required memory capacity.

Challenges in computational techniques
It is an advantageous feature that the PM can observe the temporal development of the EVDF.On the other hand, we often need only the EVDF solution in drift equilibrium.In conventional PM calculations, the equilibrium EVDF is obtained after physical relaxation of the EVDF under the electron acceleration by the external fields and scattering by gas molecules proceeding step by step every ∆t.This guarantees the continuity of the number of electrons.When the difference between the normalized EVDFs derived from f (v, t + ∆t) and f (v, t) is negligibly small, we may regard the EVDF as the converged solution.The normalized equilibrium EVDF satisfies the EBE, which is a balance equation under ∂/∂t = 0.The relaxation process of the EVDF can be accelerated by a scheme based on the Gauss-Seidel method [68].In contrast to that f (v, t) is kept unchanged until f (v, t + ∆t) is obtained in case physical relaxation of the EVDF is observed, the Gauss-Seidel method renews f (v, t) part by part (i.e.cell by cell) on the basis of the local balance expressed by the EBE ignoring the entire electron conservation.It is expected here that renewed value of f (v, t) is closer to the equilibrium value than before renewal.Fast propagation of the renewal result over velocity space is promoted by arranging the sequence of renewal calculations to be from the upstream cells to the downstream cells in velocity space.This numerical relaxation process no longer has physical meaning.However, the number of iterations for convergence to the equilibrium solution can be reduced down to the order of magnitude of 1/1000 in a drastic case [68].
A challenge to a long-term relaxation process was demonstrated with a PM in a matrix form [38].When the EVDF at a time t is represented in a vector form f (t) with the elements corresponding to the number of electrons in each cell similarly to Equations ( 17) and ( 18), f (t + ∆t) can be obtained by applying a propagator represented in a matrix form M as f (t + ∆t) = M f (t) in the same way as Equation ( 26).Here, we may calculate powers of M as by repeating the squaring operation for P i+1 = P i × P i .After this preparation, we obtain f (t + 2 n ∆t) = P n f (t).The relaxation time 2 n ∆t increases exponentially with n while computational elapse time increases linearly with n.The calculation for a slow relaxation of the EVDF in a ramp model gas [92] was demonstrated using 18 000 cells defined with 1000 divisions for ϵ ∈ [0, ϵ max ] and 18 divisions for θ ∈ [0, π] (thus, vector f (t) has 18 000 elements), ∆t = 1 ps, and n = 30 (2 30 > 10 9 , thus until t > 1 ms) [38].
Convergence of mean electron energy ⟨ϵ⟩ and the EVDF were observed midway by n = 20.This approach enables us to observe a long-term relaxation in a logarithmic time scale, although it requires a large matrix size as 18 000 × 18 000.
Another recent improvement of the PM calculation is for calculations under low E/N conditions and rf E fields [82].The polar-ϵ configuration has a tendency that the cells around the origin become more coarse than in the polar-v configuration.This feature is negligible at high E/N values where the electron energy loss by inelastic collisions and the electron energy gain under high E are dominant in the formation of EVDF satisfying the energy balance.However, at low E/N values, which appear not only under dc E fields but also in rf E fields as low E periods, the coarseness of the low-energy cells leads to an overestimation of the acceleration of low-energy electrons and causes a less precision in some electron transport coefficients.This is because uniform electron distribution within a cell is used to be assumed in conventional PM calculations.It was demonstrated that the overestimation is suppressed by blending central and upwind differences with relevantly chosen weights.
Owing to the progress in both algorithm and computational facility, the PM is nowadays recognized as an available and possible self-consistent EVDF solver.With already established techniques, the PM can be used not only to solve the EVDF in typical observation models but also for preparation of look-up tables of a set of electron transport coefficients, mean electron energy ⟨ϵ⟩, ionization coefficient ν ion , drift velocity W , and diffusion coefficient D under uniform dc E, for the use in fluid simulations.
Further enhancement of the calculation speed and memory capacity would arrow us consideration for more practical effects of electron-molecule interactions, higher-order multi-dimensional cell configurations, and more complicated geometries of the plasma reactors.Some of the basic processes not dealt with sufficiently would be anisotropic scattering, super-elastic collisions, Coulomb collisions between electrons, etc.For the cell configuration, realization of PM calculation in dimensions higher than three would have to wait for more empowerment of the computers; e.g., self-consistent PM calculation in a cylindrical symmetry requires five-variable phase space as (r, z, v, θ v , ϕ), which requires a tremendously huge memory capacity even with a coarse resolution for each variable.Nonetheless, for frequently assumed geometries such as one-dimensional real space between parallel-plane electrodes, the PM calculations done for three-variable phase space would withstand requirements on practical conditions by its flexibility for linkage with, for example, Poisson's equation and boundary conditions.

Summary and conclusions
The methods that have been described in this review have in common that they describe the kinetics of electrons, subject to the effect of electric and magnetic fields and collisions with atoms and molecules in a low-temperature plasma, using a propagatorbased description.Propagators are operators that shuffle the probability distribution of the electronic fluid in the phase space, moving it from one position to another in this generalized space: as many mathematical objects of a certain importance can have a description apparently very different from another depending on of the formal language on which the initial setting of the algorithm is based.This means that in some cases, such as in the writing of a so-called traditional MC program [27], the role of propagators may be much less evident at a superficial level.In some other cases [77], propagators enter not only into the settings of the algorithm but right from the beginning of the language that is used to describe it.A hybrid stochastic-deterministic approach, such as the MCF [37,61], which effectively uses the propagators twice, first as "micro-propagators" for calculations of transition probabilities, and then as Markov matrices to extend the simulation to the characteristic times of energy diffusion, is an even more particular case.
Emphasizing the common roots of these approaches and others that are comparable, does not erase the fact that in the concrete study of a certain type of ionized gas rather than another a method in particular may be much superior to others: elements that contribute to this choice are the relationships between the various characteristic times of the plasma, the type of collision processes, the dimensionality of the system, the presence and complexity of boundary conditions, even practical requirements such as the reduction of requested memory, speed and computer time, the physical transparency of the method and ease of code maintenance.However, the ability to understand the conceptual bases and common aspects from a mathematical point of view of the calculation methods available for the study of the electron kinetics of low-temperature plasmas can represent a very useful element of knowledge, both for researchers who have to decide the best strategy for the representation of a system and for instructors who teach these methods and arouse interest in them.
In all cases, the fundamental object in question is a time-dependent relationship between two states in phase space and which can be described physically as a function or as an operator acting on a function, and numerically as a matrix.Hence, its treatment requires tools from algebra and calculus.It is possible that the current approach to the simulation of plasma-based application systems, which is a modular approach, may hide the conceptual interest of these calculation tools: they constitute the most physical part of simulation programs and their connection with other areas of physics and mathematics are both aspects that can spark new interest in these tools.The mathematics of propagators allows us to solve the problem of the kinetics of plasma electrons with great efficiency, and this makes these instruments useful tools, because the kinetics of electrons in plasma cannot be neglected in reliably calculating chemical reaction rates and transport quantities.
Understanding the mathematical concept underlying a group of models is a stimulating exercise in itself and also has a heuristic value, because it produces future research questions: • What is the most efficient presentation of the propagation concept depending on the conditions of plasma?
• How can the concept of propagation be extended to nonlinear systems such as those that emerge when collisions between charged particles or gas heating are taken into account?
• Can the mathematical grids that are employed at different times of the different approaches be replaced with spectral descriptions?
• Could propagators, which represents input-output relations, efficiently be computed using machine learning and artificial intelligence?
• Could the similarity between the propagators used in plasma and those used in other areas of physics help to develop faster or less memory-intensive computational methods?
More recent applications of low-temperature ionized gases from biology to space science to materials science will always produce new possibilities for numerical descriptions of electron transport that have new strengths even as they have their specific weaknesses.The authors hope that this review work can stimulate research activity that captures new insights and develops new languages.

Figure 2 :
Figure 2: The MCF schematics.(a) discretization of the velocity space into cells, (b) MC simulations of electrons for calculation of transition probabilities between velocity space cells, (c) solution of Master Equation for deterministic evolution of the electron distribution function.
distribute.The number of electrons belonging to a cell is stored in an element of a matrix or an array corresponding to the cell.The probability of intercellular electron transition from a source cell to destination cells in a short time step ∆t is represented by operators called propagator or the Green function.The transition may be caused by the velocity change under electric and magnetic fields and scattering at collisions with gas molecules.The PM does not use random numbers, and its calculation is completely deterministic.The stochastic processes are condsidered on the bases of their expected values.Temporal development of the EVDF is observed by applying the acceleration and collision propagators to the EVDF every ∆t.The PM has a flexibility in treatment of the propagators as seen in practical examples introduced in the following sections.The step-by-step time development of the EVDF can be modified in some specific models customized to derive only equilibrium solutions, and the PM can deal with some simple boundary conditions such as electron reflection and secondary electron emission at electrode surfaces.

For
the latter point, Lagrangian and Eulerian expressions are considered for the treatment of electron motion in the acceleration or flight.In the Lagrangian expression, the source and destination cells are related by ballistic motion of electrons starting from the source cell (Figures3a and 3c), Think that the electrons in a source cell undergo a collective free flight for ∆t.They would appear as if they are a moving cell, which is called a Lagrangian cell.The destination cells are chosen as those located at the position where the Lagrangian cell reaches after ∆t, and the ratio of the electron redistribution that the destination cell accepts from the source cell is calculated as the ratio of the area overlapping between the Lagrangian cell and each destination cell.On the other hand, in the Eulerian expression, source and destination cells are assumed to contact each other and to share the cell boundary between them.The number of electrons moving out of the source cell to the destination cell, n e,out , is evaluated by integrating the electron flux passing through the cell boundary during ∆t (Figs.3b and 3d).Its amount is evaluated as n e,out = S cell (n e,cell /V cell )(eE/m)∆t, where S cell [L 2 T −2 ] is the area of the cell boundary projected to a plane perpendicular to E, n e,cell is the number of electron in the cell, V cell [L 3 T −3 ] is the volume of the cell, and (eE/m) [LT −2 ] is the electron acceleration.∆t must be short enough to avoid n e,out > n e,cell for all cells.

Figure 3 :
Figure 3: Cells defined in two-variable velocity space (v ∥ , v ⊥ ) or (v, θ), and treatment of the intercellular electron transition by acceleration: a, Lagrangian treatment in Cartesian-v cell configuration; b, Eulerian treatment in Cartesian-v cell configuration; c, Lagrangian treatment in polar-v cell configuration; and d, Eulerian treatment in polar-v cell configuration.The thick boundaries indicate the source cells from which electrons flow out downstream (rightward) by acceleration.The red hatched cells are the Lagrangian cells.

Figure 4 :
Figure 4: Treatment of electron collisions and scattering: a, elastic collision, electron redistribution to the cells of ϵ = ϵ ′ ; b, excitation collision, electron redistribution to the lower-energy cells of ϵ = ϵ ′ − ϵ exc ; c, ionization collision, electron doubling and redistribution to the lower-energy cells of ϵ ≤ ϵ ′ − ϵ ion ; and d, electron attachment, electron disappearance from velocity space.The cells with thick boundaries and red crosses are the source cells, from which electrons flow out.ϵ ′ is the energy of the incident electron.ϵ exc and ϵ ion are the excitation and ionization energies, respectively.
sin θ are components of v parallel and perpendicular to the direction of −E, respectively, and v = |v|.
is the weight of the kth-order Hermite function, G(t) = ⟨z⟩ is the center of mass of the whole electron swarm, σ(t) is the standard deviation of electron position z, and Z = z/[ √ 2σ(t)] is the dimension-less z position.H k (Z) satisfy an orthogonality 81,87,88],87,88]l variation of G(t) is the centerof-mass drift velocity W r (t) = dG(t)/dt of the electron swarm[68,80,81,87,88].m 2 (v, t) gives the variance of the electron swarm as [ v m 2 (v, t)dv]/[ v m 0 (v, t)dv] − [G(v, t)] 2 .Its temporal variation gives the longitudinal diffusion coefficient as DL (t) = 1 2 (d/dt){ v m 2 (v, t)dv/ v m 0 (v, t)dv − [G(t)] 2 } [68,81,87,88].For the moment calculation up to the nth order, (n + 1) sets of cells are prepared, and initial values ofm k (v, t) (0 ≤ k ≤ n) arestored in the cells.Their temporal variations for ∆t are calculated by applying the collision and acceleration propagators to m k (v, t).The propagators are common for m k (v, t) irrespective of k.In addition, amounts of the kth-order moments kv z m k−1 (v, t)∆t evaluated by the drift term are added to m k (v, t) to obtain m k (v, t + ∆t).The instantaneous amount of the kth order moment of the whole electron swarm is given as m k (t) = v m k (v, t)dv.m k (t), which are values in laboratory system, are converted to the values in center-of-mass system around G(t) as m ′ k are further converted to to the dimension-less valueM k (t) reduced by σ(t) = m 2 (t)/m 0 (t) − [G(t)] 2 as M k (t) = m ′