Temperatureinduced negative differential conductivity
Abstract
We report on the existence of temperature induced negative differential conductivity (NDC) for electrons in gaseous nitrogen. The important role of superelastic rotational collisional processes in this phenomenon is highlighted. A model crosssection set, utilised to ensure an accurate treatment of superelastic processes and achieve thermal equilibrium is detailed, and used to illustrate the role of deexcitation processes in NDC. The criterion of Robson (Robson, 1984) for predicting the occurance of NDC using only knowledge of the collisional crosssections is utilised for both the model system and N We also report on the impact of anisotropy in the very low threshold scattering channels on the transport coefficients, examine the FrostPhelps finite difference collision operator for the inelastic channel, in particular its neglect of recoil, and assess other assumptions utilised in existing Boltzmann equation solvers. We discuss the numerical challenges associated with low reduced electric field calculations, and detail an alternative representation of the elastic and inelastic collision operators used in Boltzmann equation solutions that enforce conservation of number density. Finally, new experimental measurements of the drift velocity and the Townsend ionisation coefficient for an electron swarm in N are reported from a pulsed Townsend experiment. The selfconsistency of the utilised crosssections is also briefly assessed against these results.
I Introduction
As a swarm of electrons drift and diffuse through a background medium, driven out of equilibrium by an externally applied electric field, there can exist a region where the drift velocity of the electrons decreases with increasing electric field strength. This phenomenon, known as negative differential conductivity (NDC), has been comprehensively studied, both experimentally and theoretically (Petrović et al., 1984; Dyatko et al., 2014). In both plasma and swarm physics NDC is present in gases used for dosimetry and particle detectors (Christophorou et al., 1979; Mathieson and El Hakeem, 1979; AlDargazelli et al., 1981), and has implications on the operating ranges of gas lasers with NDCinduced electric current oscillations in electronbeamsustained discharge switches (Lopantseva et al., 1979; Christophorou, 1981). For fundamental physics, NDC has played a role in evaluating complete and accurate scattering crosssection sets (Petrović et al., 2007; Nakamura, 1995). Argon, for example, was considered to be a candidate for NDC in a pure gas, but this was later shown to be due to the presence of molecular impurities in Ar samples (Long Jr et al., 1976), and in gaseous mercury it has recently been shown that NDC occurs due to the presence of dimers (Mirić et al., 2017). Swarms of electrons can also induce NDC in liquids (White and Robson, 2009), plasmas and semiconductors (Chiflikian, 1995; Chiflikyan, 2000; Dyatko, 2007), and NDC has been shown to be induced by positron swarms in argon (Šuvakov et al., 2008). In addition to NDC for positrons in argon, the same phenomenon has been observed in water vapour (Banković et al., 2012a), molecular hydrogen (Banković et al., 2012b) and CF (Banković et al., 2014). As such, modelling systems to predict regions of NDC and the conditions of electroninduced NDC is of particular interest (Robson, 1984; Petrović et al., 1984).
In model systems, the work of Petrović et al. (Petrović et al., 1984) and Vrhovac and Petrović (Vrhovac and Petrović, 1996) detail different systems involving elastic, inelastic and ionisation crosssections that, under various conditions, either enhance or eliminate the occurrence of an NDC region. Further to the electron and position experimental studies, in real systems, the existence of NDC has been observed experimentally in gaseous systems of N (Pack and Phelps, 1961a), CH, CF (Hunter et al., 1985) and Hg (Mirić et al., 2017), and predicted theoretically due to electronelectron interactions ((Chiflikyan, 2000)) in plasmas of Xe (Aleksandrov et al., 1996; Donko and Dyatko, 2016), to name but a few. NDC in gas mixtures also has an extensive history, observed in mixtures with helium, argon, N and CH (Chiflikian, 1995; Hunter et al., 1985; Dyatko et al., 2010, 2014; Stano et al., 2011). In strongly attaching gases, NDC has also been shown to be induced through a combination of attachment heating and inelastic cooling (Mirić et al., 2016).
Throughout this broad body of work, the sources of NDC have been discussed in detail by the various authors and a number of criteria have been proposed. Some of the conditions under which NDC can occur include the presence of inelastic collision channels, favoured particularly by a decreasing inelastic crosssection, the presence of a RamsauerTownsend minimum in the elastic momentumtransfer crosssection, or a rapidly increasing elastic momentumtransfer crosssection. However these are not necessary and sufficient conditions. The validity of these early criteria on the understanding of NDC was discussed in detail in Petrović et al. (Petrović et al., 1984) who address the analyses of Kleban and Davis (Kleban and Davis, 1977, 1978), Long and coworkers (Long Jr et al., 1976), and Lopantseva and coworkers (Lopantseva et al., 1979), in particular, and in Vrhovac and Petrović (Vrhovac and Petrović, 1996) where consistency with the Shizgal (Shizgal, 1990) criterion is discussed. Of particular note are the simulations of Petrović et al. (Petrović et al., 1984) for N that confirm that the presence of inelastic processes, other than rotational excitations, decreases the range of the NDC region, indicating that the rotational collisions are responsible for the presence of NDC, where the elastic crosssection is relatively isotropic. This is explored further in this work, where the effect of temperature on the presence or absence of NDC is explored. The thermally induced NDC region in N below room temperature is detailed. Here, the contribution of inelastic groundstate and excitedstate collisional processes to the nett energy transfer to and from the electron swarm were shown to be responsible for the extent of an NDC region.
The various criteria for the presence of NDC, detailed in these early works, as noted above, have been discussed in detail by the respective authors. Of particular interest here is the criterion proposed by Robson (Robson, 1984), where momentumtransfer theory was used to derive an expression based on the rate of change of the ratio of inelastic to elastic energy transfer with energy. This criterion allows prediction of NDC using only a knowledge of the collisional crosssections, the accuracy of which is highlighted using both a simple model system and for N. The thermallyinduced reduction and deactivation of NDC is explored further using Robson’s criterion.
A model collisional system is used throughout this study to simplify discussions around NDC and its temperature dependence. The model system is also used to verify both our solution method, and explicitly test the inclusion of temperature dependent inelastic collisions. This system also facilitates further discussion around some of the assumptions sometimes involved in swarm modelling. The explicit effect of anisotropy in lowthreshold inelastic processes is assessed, as well as the commonly employed twoterm approximation (White et al., 2003), as is the neglect of superelastic collisions in higherorder collisional terms. Also assessed is the neglect of the recoil of the neutral particle during inelastic collisions in the excitation collision operator, as utilised in many Boltzmann equation solutions. Truncation of the mass ratio expansion at zeroth order for inelastic collisions is compared with the exact collision description of Monte Carlo simulations, for very lowenergy threshold processes like rotational excitations, where the energy exchange is much closer to that for elastic collisions. For application of this discussion to real gases, these assumptions are also assessed for electron swarms in N. For the low energy regime of interest in this study the discrete rotational collisions in N are treated using the FrostPhelps collision operator, although in the work of Ridenti et al. (Ridenti et al., 2015) the ChapmanCowling extension to the continuous approximation to rotations was developed to bridge the continuous energy loss regime applicable at high fields to the discrete collision description for use at low energies.
The particular numerical techniques and challenges associated with modelling lowthreshold processes, in low electric field conditions, are discussed for both our model system and in N. For a robust solution method to reproduce various experimental conditions, standard collisional benchmarks are employed along with selfconsistency checks of numerical number density conservation, energy and momentum balance. Under particular conditions we are able to provoke nonphysical solutions. The effect of these problematic solutions is assessed and addressed using alternative representations of the standard collision operators, and compared with the results of an independent Monte Carlo code.
This work is arranged by first detailing our kinetic theory, our multiterm solution and particular numerical approaches, in section II and Appendix A. The Monte Carlo code used for comparison with our results is briefly outlined in section II.5, while the pulsed Townsend apparatus, used for measurement of the drift velocity and the Townsend ionisation coefficient in N, is also briefly detailed in section II.6. In section III we consider two systems that exhibit NDC—a simple model system in section III.1 and N in section III.2. For both systems we present calculations of transport coefficients at various temperatures, and discuss the thermal activation of NDC and the necessity of deexcitation, or superelastic, collisions for its occurrence. The prediction of NDC using Robson’s (Robson, 1984) criterion from knowledge of the energy transfer rates is presented for both systems, in sections III.1.1 and III.2.4. We also consider various physical and numerical elements including: (i) the effect of the neglect of recoil in the FrostPhelps differential finite difference inelastic collision operator, (ii) the twoterm approximation on the Townsend ionisation coefficient, (iii) the effect of anisotropy in lowthreshold inelastic collision channels of our model crosssection, and (iv) the neglect of higher order deexcitations sometimes applied in the solution of similar problems, in sections III.1.3 and III.2.3. In section III.2.2 and Appendix B, we also present new experimental measurements of the drift velocity and the Townsend ionisation coefficient for electron swarms in N to compare against our calculations.
Ii Theory
ii.1 The Boltzmann equation and a multiterm solution framework
The transport of a swarm of charged particles through a gaseous medium is described by the particle phase space distribution function , representing the distribution of electrons with position , velocity and time , that is the solution of the linear Boltzmann equation,
(1) 
where is the mass of the swarm particle, is the elementary charge, and is the externally applied electric field. The linear collision operator describes binary collisions between the swarm particles and the background medium, and accounts for elastic and inelastic collisions, and particle nonconserving loss (attachment) and gain (electron impact ionisation) collisions.
For a solution to the Boltzmann equation, the angular dependence of the phase space distribution function is expanded in terms of spherical harmonics, to give:
(2) 
In practice, the index must be truncated at some upper value , incremented until some convergence criterion on the distribution function, or its velocity moments, is met. For this work we do not restrict the truncation at , as is commonly done for the ‘twoterm approximation’, sometimes leading to an inadequate representation of the anisotropic parts of the distribution function and incorrect transport coefficients, see the review (White et al., 2003). Substitution of the expansion into equation (1) leads to a system of coupled equations for (Robson et al., 2002).
For the plane parallel geometry representing our experimental conditions, the preferred direction is taken to be perpendicular to the electrodes and the spatial gradients are along the axis so that and , and the index is restricted to by symmetry. Equation (1), with substitution of the expansion (2) and recast in energyspace using for in eV, becomes:
(3) 
where is the Legendre decomposition of the collision operator, detailed in the next section, and is the electric field defined parallel to the axis.
ii.2 Collision operators
ii.2.1 Elastic collisions
For electron swarms in atomic and molecular gases the small mass ratio is utilised so that the Davydov operator for elastic collisions may be used (Davydov, 1935; White and Robson, 2011; Pidduck, 1936), and is given by:
where is the momentumtransfer collision frequency, is the th partial collision frequency, in , is Boltzmann’s constant, and is the temperature of the neutral background gas.
ii.2.2 Inelastic collisions
The Frost and Phelps Legendredecomposed collision operator (Frost and Phelps, 1962) is employed here to describe the effect of inelastic particleconserving collisions on the spatiallyindependent velocity distribution function. The anisotropic form of the collision operator was detailed by Makabe and White (Makabe and White, 2015), Phelps and Pitchford (Phelps and Pitchford, 1985) and earlier in Reid (Reid, 1979) (in the second and third terms on the right hand side of the equation following equation (3)), and is here extended to include deexcitation, or superelastic, collisions using detailed balance (Hochstim, 1969). Expressed in terms of initial and final internal states and of the neutral particle, where , particles with energy above the threshold are available for excitations from . Below the threshold, for nonzero temperatures, the background neutral particles may be in an excited state and are available to undergo superelastic collisions from , where the energy loss, taken to be the threshold, is gained by the incoming electron in a deexcitation process and lost by the neutral particle.
The partial crosssections are the coefficients of a Legendre polynomial () expansion of the differential crosssections, for the scattering angle , defined by . The deexcitation crosssection is expressed in terms of the excitation crosssection using the microscopic reversibility relation where and are the degeneracy of the th and th states. After converting to collision frequencies through in energy space, where is the neutral number density, the isotropic and anisotropic components of the inelastic collision operator are given by:
where and are the density of neutral particles in the initial and final states and , respectively, is the th partial collision frequency, related to the momentumtransfer crosssection through for inelastic collisions. The number density of the neutral particles in the state with energy , are calculated using standard Boltzmann statistics: where the partition function sums over all possible internal states and is given by .
ii.2.3 Ionising collisions
The electronimpact ionisation operator utilised here, in Legendredecomposed form, is given by:
where is the ionisation collision frequency, and is the energypartitioning function (Ness and Robson, 1986).
ii.3 Solution technique
Timeofflight:
When the number density varies slowly in space, away from boundaries and under the influence of a uniform electric field, hydrodynamic conditions prevail and the spacetime dependence of the distribution function can be projected onto the number density so that , where is the rank of the tensor.
To account for particle nonconserving processes, the density gradient expansion to second order is required. In planeparallel geometry this is given by:
where the and superscripts on the distribution are defined parallel and transverse to the electric field, respectively. For weak gradients, a density gradient expansion of the phasespace distribution function may be taken and the resulting diffusion equation (Robson et al., 2011):
is used to analyse experimental parameters. The timeofflight coefficients of the density gradient expansion are found from the solution to the hierarchy,
(6)  
where and the are given by,
The equation for the first level transverse distribution function takes a different form:
Each of the expansion coefficients satisfy the normalisation condition . The coefficients , using the density gradient expansion, are given by:
(7)  
where is the particle nonconserving, or reactive, collision operator. In the presence of nonconservative collisions these coefficients involve an integration over the density gradient expansion coefficients , so that these expressions become nonlinear.
Steadystate Townsend:
In modelling the steadystate Townsend experiment (Huxley, 1940; Townsend, 1915; Huxley and Crompton, 1974), for the steadystate timeasymptotic solution, far from the source where no memory of the initial source distribution remains, the spatially varying distribution function takes the form of a sum of exponentials , where and are separation constants (Robson et al., 2011). Substitution into equation (3), leads to the generalised eigenvalue equation:
(8) 
where the coefficients are defined in equation (3) and the eigenvalues and are related through the dispersion relation .
For the steadystate Townsend experiment, the spatial eigenvalues are found by setting to zero, which implies we are looking at the temporally asymptotic regime where . The spatial eigenvalues are assumed to form a discrete set where the lowest magnitude nonzero eigenvalue represents the reduced macroscopic Townsend ionisation coefficient .
Utility of the generalised eigenvalue method:
The generalised eigenvalue equation (8) may be solved directly as an eigenvalue problem, the details of which are omitted here but are described in Boyle (Boyle, 2015). With this technique, the lowest remaining temporal eigenvalue of a timeofflight simulation represents the rate of nonconservative processes in the longtime limit, equivalent to in equation (7), and the corresponding eigenfunction represents the electron energy distribution function. This provides an alternative method for calculating the nett rate coefficient of nonconservative collisional processes for a timeofflight simulation that is typically calculated as the integral of the nonconservative collision frequencies with the distribution function, defined below in equations (9)–(10).
Numerics and benchmarking:
The particular numerical methods employed in the solution of equations (6) and (8) are detailed in reference (Boyle, 2015). Here a nonuniform energy grid is employed, that is dense at lower energies to capture the variations near the lowenergy thresholds of the inelastic processes considered throughout this work. The other explicit changes to the methods of (Boyle, 2015) used here are described in Appendix A. Systematic benchmarking of the theory and numerical solution has been performed.
ii.4 Transport coefficients
Knowledge of the full phasespace distribution function allows for the calculation of all macroscopic quantities describing the electron swarm. The distribution function , that is the solution to equation (6), allows the calculation of the timeofflight transport coefficients. The coefficients used here include the mean energy , flux and bulk drift velocities, the nett rate coefficient summed over all reactive collision frequencies , and the bulk longitudinal diffusion coefficient , and are defined as:
(9)  
(10)  
In the absence of nonconservative collisions, the bulk transport coefficients reduce to the flux coefficients, and , represented by the first term in each of the bulk coefficient definitions.
An important selfconsistency/accuracy check for an accurate solution are the rates of energy and momentum exchange, where the gain from the advective terms (the external electric field and time rate of change components) and loss due to collisions must be balanced. Calculation of energy and momentumtransfer rates due to individual crosssections allows assessment of the contribution of not only each collision type, but separation into inelastic and superelastic channels. For the timeofflight experiment, the power exchange due to the advective terms is given by:
while the power exchange due to each collisional process is given by:
For the steadystate Townsend configuration, denoted by the subscript , the mean energy and drift velocity are given by (Robson, 1991; Dujko et al., 2008a):
where is the spatially averaged mean energy, is the Townsend ionisation coefficient, is the gradient energy parameter (White et al., 1995), and is the flux longitudinal diffusion coefficient given by the first term in the bulk coefficient definition. In the absence of nonconservative collisions, the coefficients reduce to the flux coefficients. The Townsend ionisation coefficient, calculated directly from the solution to equation (8), can also be related to the bulk hydrodynamic timeofflight coefficients, when spatial gradients are weak, through the nett reaction rate (Robson, 1991; Dujko et al., 2008a):
(11)  
(12) 
where the summation is over all of the reactive processes , like ionisation and attachment, and the expression for the Townsend ionisation coefficient to second order becomes . The experimentally measured Townsend ionisation coefficient is related to the calculated coefficient through .
ii.5 Monte Carlo technique
We have implemented a standard swarm MonteCarlo sampling code. The code uses the nullcollision method (Skullerud, 1968) along with temperature included via appropriate modifications to the total crosssection and resolution of collisions (Ristivojevic and Petrović, 2012). Measurements are made through the ‘box sampling’ style (Dujko et al., 2008b), where an integral over the quantities to be measured is performed between each collision and binned into time bins. Hence a timespecific measurement refers to an average of that quantity during the time bin. To ensure we have considered a large enough simulation time to have reached steadystate, we consider a sufficiently fine time grid to allow a fit of the quantities to the empirical form of: , where is the steady state value for quantity and is the final time of the simulation. This definition allows us to give the meaning of a deviation from steadystate at the halfway point of the simulation. The condition, is enforced, and we then average over the latter half of the simulation to build up the statistics for the MonteCarlo results.
We estimate the error in these results by the standard error of the averaged simulations at different times. We have also ensured that the autocorrelation between consecutive points is minimal.
The MonteCarlo code has been tested against many benchmarks including pure elastic models of hard sphere and Maxwell models (Boyle, 2015), argon measurements (Boyle et al., 2015), the inelastic and anisotropic models of Reid (Reid, 1979; Boyle, 2015), the ionisation models of Ness and Robson (Ness and Robson, 1986; Boyle, 2015), and inclusion of a static structure factor (Tattersall et al., 2015).
As part of the tests to be performed in section III, we require a different temperature for the elastic and inelastic processes. We have implemented this by considering a mixed system of two species. The first species possess only an elastic process, with a gas temperature given by the elastic temperature. The second species possesses only an inelastic process, with the ground and excited populations given by the inelastic temperature. When the elastic and inelastic temperatures coincide, this is equivalent to a simulation of a single species with both processes.
ii.6 Experimental technique
The fully automated pulsed Townsend experiment, used to measure the drift velocity and the Townsend ionisation coefficient for electrons in gaseous , has been described in detail previously (HernándezÁvila et al., 2002; Basurto et al., 2013; de Urquijo et al., 2014) and so the experiment is only briefly summarised here. The total displacement current of the electrons, and their ionic products, that drift through the parallel plate capacitor under a homogeneous electric field is measured and separated into a fast component due to the electrons, and a second part due to the slower ions.
The initial swarm of electrons is generated from the cathode by an incident 3 ns duration, UV (355 nm) laser pulse, and the electrons and ions formed by reactions with the neutrals drift to their respective electrodes under the action of a highly homogeneous electric field , produced by a very stable voltage in the range 0.25 kV, according to the densitynormalised field selected and the density of the gas in the discharge vessel. The electrons and ions drift through the capacitor with a fixed drift distance of 3.1 cm (0.025 mm), between an aluminium cathode and a nonmagnetic stainless steel anode, each 12 cm in diameter. The electrons are detected with a lownoise, 40 MHz amplifier with a transimpedance of V/A. The measurements presented here were performed over the temperature range 293300 K, measured with a precision of 0.5 K, and with a pressure range of 0.530 Torr, as monitored with an absolute pressure capacitance transducer with 0.15% uncertainty. The commercial grade sample of N used here from Praxair had a stated purity of 99.995%. The overall uncertainty in the measurements of the electron drift velocity was 2.2%, and 7.4 to 9.4% for the Townsend ionisation coefficient.
Iii Results and Discussion
iii.1 NDC — a model crosssection study
A limitation of the existing collision benchmark models we use is in testing/verifying the inclusion of superelastic processes in the inelastic channel. In the absence of superelastic processes, thermal temperatures cannot be achieved, so a simple model system, verified by an independent Monte Carlo method, allows us to confirm that our solution methods are well representing physical processes.
In the pursuit of clear criteria for the existence or prediction of NDC, a number of model crosssections have been proposed (see for example (Petrović et al., 1984)). Many of these models could be adapted to account for the inclusion of superelastic processes. The model considered in this work, however, was chosen to illustrate the damping effect of superelastic populations on NDC at room temperature, similar to the behaviour of electrons in molecular nitrogen. For collisions with neutrals with a mass at 0 K, 77 K, and 293 K the transport coefficients have been calculated for the model elastic and excitation crosssections (in atomic units):
(13) 
where and .
Figure 1 shows the drift velocity and mean energy calculated using the multiterm Boltzmann equation solution and the independent Monte Carlo code, as a function of the reduced electric field in units of the Townsend (1 Td = ). The agreement between the Monte Carlo and Boltzmann solutions is better than 2% for the drift velocities and the mean energies at 0 K, and with less than a 4% variation in the drift velocity and a 2% variation in the mean energy at 77 K. However, this increases to 6% for both the drift velocity and mean energy at 293 K.
To address the discrepancy between the transport coefficients calculated using our Boltzmann and Monte Carlo solutions, we consider the neglect of recoil in the inelastic channel in our solution (and similar solutions of the Boltzmann equation, for example the recent work of Ridenti et al. (Ridenti et al., 2015) in the continuous energy loss approximation). Unlike elastic collisions, that are represented to first order in the mass ratio to take into account the thermal motion and recoil of the neutral particle during an elastic collision, recoil of the neutral particle during inelastic collisions is neglected in most of the existing Boltzmann equation solutions, and in our solution. This assumption has been considered previously in White et al. (White et al., 2002), using the integral form of the inelastic collision operator that does not restrict collisional representation to zerothorder. The transport coefficients, calculated for electron impact on H using a multiterm solution over the range 0.1–10 Td, differed by less than 0.1% between no recoil inelastic collisions and the converged collision description. In their work, the lowest excitation channel in H is the rotational excitation with a threshold of 44 meV, while for our model system the energy loss threshold is 2 meV, much closer to the first order mass ratio meV for the model system.
To include the thermal motion of the neutrals during inelastic collisions in the FrostPhelps differential finite difference collision operator, requires an extension that is outside the scope of this study. However, we still desire quantification of the effect on the transport coefficients. In Monte Carlo simulations the collisions are treated exactly, so the Monte Carlo technique described in section II.5 was used to assess the effect of recoil in the lowthreshold channel of interest here, as shown in figure 1. The effect of truncation of the mass ratio for inelastic collisions on our calculations is most prevalent at reduced electric fields between 0.1 Td and 10 Td, where the nett energy transfer due to elastic collision is increasing relative to the energy transfer due to inelastic collisions, as shown in figure 2. Here, the difference between the complete and approximate collision descriptions in the Monte Carlo calculations is greatest at 293 K, where the drift velocity and mean energy both differ by up to 7%. At 77 K, the differences are up to 4% and 3% between the drift velocity and mean energy, respectively, and at 0 K the drift velocity and mean energy differ by 1.7% and 1.3%, respectively, between the two collision representations.
When recoil of the neutral particle during inelastic collisions in our Monte Carlo simulation is neglected by artificially increasing the neutral mass for inelastic collisions only, to replicate the differential finite difference form of the inelastic collision operator utilised here, the difference between our Monte Carlo and Boltzmann calculations reduces to 0.9% and 1.3% for the drift velocity and mean energy, respectively, at 0 K, 1.4% and 0.8% between the drift velocity and mean energy at 77 K, and generally below 2.8% (increasing to 4.7% at low fields with statistical noise) and 1.1% between the drift velocity and mean energy at 293 K, respectively.
The presence of NDC is anticorrelated with the presence of superelastic collisions, highlighting the damping effect of the deexcitation process on the presence of NDC. When properly included through detailed balance for inelastic collisions, a smaller ratio of neutrals in the groundstate, caused by an increasing temperature, increases the mean energy of the swarm and decreases the drift velocity. The energy transfer profiles given in figure 2 show the lower nett inelastic energy transfer rate with increasing temperature, increasing the mean energy which samples higherenergy regions of the elastic crosssection, resulting in a reduced average velocity of the swarm. NDC ceases when the collisional energy transfer is dominated by the elastic process. At higher temperatures, the increased fraction of neutrals in excitedstate populations reduces the nett power transfer due to inelastic collisions, as shown by the superelastic contribution to the energy transfer in figure 2. As a direct result of the superelastic population, the range of NDC is reduced and the transition to elasticcollision dominated energy transfer occurs at lower reduced electric fields for increasing temperatures.
iii.1.1 Criterion for NDC
Using momentumtransfer theory, Robson’s (Robson, 1984) criterion for the presence of NDC uses the energy variation of the ratio of the elastic to inelastic energy transfer. For a monotonically increasing elastic collision frequency and open inelastic channels, the mean energy increases with increasing reduced electric field, slowly as the inelastic collisions take energy from system, so that the drift velocity increases with field, as illustrated in figure 1. As the inelastic collisions become less important relative to the elastic collisions, the mean energy of the swarm increases at a greater rate, to sample the higher elastic collision frequency, causing the drift velocity to begin to decrease with increasing field. The criterion derived for the appearance of NDC by Robson (Robson, 1984), at a given mean energy , is given by where represents the ratio of the total inelastic to elastic energy transfer and is given by:
(14) 
where the total inelastic energy transfer is taken as the sum over all inelastic channels with associated thresholds .
We note that the NDC criteria of Robson (Robson, 1984) and Petrović et al. (Petrović et al., 1984) differ due to a different expression for the energy balance, where the latter omit energy transfer due to elastic collisions, although this is sufficient for the systems considered in that work.
The criterion proposed by Robson (Robson, 1984) is a very good predictor for NDC for the model crosssection considered in this study, as shown in figure 3. The NDC region for 0 K and 77 K ceases when the energy transfer rate due to elastic collisions is greater than the nett energy transfer rate due to the inelastic process, as predicted. For the calculation of at 293 K, for reduced electric fields between 0.1 Td and 0.2 Td, the derivative is and the presence of NDC is only weakly predicted, but does not occur in our calculations.
iii.1.2 Temperature dependence and detailed balance
To illustrate the physical dependence on the inclusion of superelastic collisions, in this subsection we consider two unphysical modifications to our model system that address the effect of temperature from each of the scattering channels separately. The two models considered incorporate different temperatures for the background gas through elastic collisions and excited state populations.
The first model includes temperature dependence of the excited state population, but considers elastic collisions with stationary neutrals, equivalent to a temperature of 0 K in the elastic collision operator, with the notation K, K in the following figures. The second model involves elastic collisions with nonstationary neutrals, at 293 K and 77 K, but inelastic collisions from ground to excitedstates only, with the notation K, K in figures 4–6.
Figure 4 displays the transport coefficients from our Boltzmann solution, with a zeroth order mass ratio representation in the inelastic collision integral, and compares them with our independent Monte Carlo solution, with inelastic collisions treated exactly, for each of these models. For the model K, K, differences of less than 6% and 5.4% were found between the drift velocity and mean energy, respectively. For the two models with elastic collisions taken to be with nonstationary neutral particles and inelastic collisions between groundstate neutrals only, we find differences of less than 2.5% and 1.4% for elastic collisions at 293 K, and 2% and 1.3% when elastic collisions are taken at 77 K, between the drift velocity and mean energy, respectively.
For properly included superelastic collisions, but when temperature is not included in the elastic channel, the mean energies approach the appropriate thermal value of with decreasing reduced electric field. At 293 K a difference of around 3.6% between the model and standard calculation is observed, generally decreasing with increasing reduced electric field strength. We find that the drift velocity calculations lie very close to, but just above, the standard model calculations at 293 K, by at most 2.2%. While those differences are not as significant as some of the others discussed in this study, these calculations do illustrate the importance of detailed balance in swarm calculations for achieving correct thermal distributions.
A more dramatic difference is observed when temperature effects are taken into account for the elastic collisions, but not in the inelastic channel. Here, the calculated drift velocity and mean energy approach the 0 K calculations due to the dominance of the inelastic channel at low reduced electric fields, as can be seen in the energy transfer profiles given in figure 5. The difference between these models and the 0 K profiles shows the explicit contribution of the temperature term in the elastic collision operator. The variation from the standard temperature treatment profiles is large, and results in an overestimate of the drift velocity below 30 Td and the presence of an NDC region that is larger than that at 77 K and absent from the 293 K calculations.
When temperature effects are included through the elastic collision operator and the inelastic groundstate density, but detailed balance is not achieved due to the neglect of superelastic collisions altogether, we observe a dramatic effect on the transport coefficients, given the large contribution of the deexcitation process to the energy transfer. Although not shown, the resulting drift velocity and mean energy profiles lie between the 0 K and 293 K results, as is expected with less energy lost in the inelastic channel than the 0 K simulation, but a greater nett energy loss in the inelastic channel than the 293 K simulation, where the deexcitation collisions contribute to energy gained by the electron swarm. These model systems highlight the necessity of detailed balance in collisional processes when modelling real gaseous systems in the low energy regime.
For these nonphysical models, which disregard thermal effects in the elastic, inelastic and superelastic channels, NDC is present in our results for the two cases where the deexcitation process is removed. Similar to the results shown above, the presence of superelastic collisions increases the mean energy of the electron swarm, increasing the drift velocity monotonically. As demonstrated by the earlier results, Robson’s criterion for the presence of NDC gives a very accurate prediction based only on a knowledge of the energy transfer rates. Illustrated in figure 6 is the rate of change with energy of the ratio of the inelastic to elastic energy transfer rate, , where the rapid decrease in this rate to below 1 at around 0.1 Td corresponds to the start of the NDC region for both of the models with no superelastic processes. Regardless of the temperature of the model system, in the absence of the deexcitation process the ratio of the energy transfers in equation (14) decreases more rapidly with energy than it would otherwise, resulting in NDC until the elastic energy transfer rate starts to dominate.
iii.1.3 Approximation effects: Anisotropy in the inelastic channel and higherorder superelastic terms
In this study we are interested in the effect of anisotropic scattering in lowthreshold processes like rotational excitations. The impact on the calculated transport coefficients are explored in this subsection, alongside our assessment of other assumptions, including the twoterm approximation and those associated with the inclusion of superelastic collisions.
The effect of anisotropic scattering in the inelastic channel has been recently investigated in Janssen et al. (Janssen et al., 2016), for an excitation scattering channel of (simplified) argon with a threshold at 11.828 eV. We note this energy threshold is much higher than the lowest inelastic thresholds of molecules like N. For the low threshold processes considered here, to quantify the effects of the anisotropic terms in the inelastic operator, given in equations (II.2.2) and II.2.2, we introduce an angular scattering component for our model excitation in order to emulate the forwardpeaked nature of rotational excitation. Using a forward scattering model where the differential inelastic crosssection is , the inelastic partial crosssections are given by:
(15) 
To test explicitly the assumption of isotropic scattering in the inelastic channel, here we modify only the inelastic momentumtransfer, leaving the elastic momentumtransfer crosssection fixed, as has been considered previously by Reid (Reid, 1979) and Phelps and Pitchford (Phelps and Pitchford, 1985), for example. Note that this does not fix the total momentumtransfer crosssection. Our calculated drift velocities for the anisotropic model combining equations (13) and (15) are shown in figure 7. For the various temperatures considered in this work, the effects of anisotropy in the inelastic channel are greatest where momentum exchange is dominated by the inelastic channel. At 0 K, this difference occurs over the range 0.02 to 0.2 Td with a variation of less than 5% in the drift velocity and less than 9% in the mean energy of the swarm (not plotted). For the 77 K and 293 K simulations, the maximum difference occurs at the lowest reduced electric fields, where the momentum exchanged during superelastic collisions increases the total momentum exchanged in the inelastic channel. At 77 K the drift velocity and mean energy differ by 11% and 6.5% respectively, and 10% and 4% at 293 K, respectively, decreasing with increasing reduced electric field for both temperatures.
The validity of using a twoterm approximation has been discussed previously (e.g., in (White et al., 2003; Phelps and Pitchford, 1985)), and we briefly consider the effects of that assumption on our model calculations here. Figure 7 shows our calculated drift velocity for isotropic and anisotropic scattering using the model crosssections of equations (13) and (15). For both the isotropic and anisotropic models, the difference between the twoterm and multiterm results is greatest at 0 K, with a difference of up to 8% for reduced electric fields below 0.1 Td in the drift velocity, and 5% in the mean energy. While at 77 K and 293 K, for both models, the differences between these transport coefficients reduces to below 0.3%.
The final assumption we wish to address is that while deexcitation is considered in the equation of the inelastic operator, it is sometimes neglected in the equations. For our isotropic model detailed in equation (13), we have removed superelastic collisions in the channels in two ways. First we consider the proportion of particles in the groundstate calculated according to the neutral temperature, as is included through the equation, but simply turn off the deexcitation channel, denoted by the notation , for example, in figure 8. The differences between the calculated drift velocity and mean energy, when compared with our standard treatment, are up to 33% and 16% at 77 K and 30% and 10% at 293 K, respectively. We also consider neglecting higherorder superelastic terms by setting all neutral particles in the ground state, with the notation , and find much smaller differences as the total number of excitation collisions remains constant, with differences of less than 9% and 1% in the drift velocity and mean energy, respectively, at 77 K, and 2% and 0.3% at 293 K. The magnitude of these differences decreases with increasing reduced electric field, influenced by the relative strength of the two crosssections, and the dominance of the elastic crosssection above 1 Td. We have also repeated these calculations using our anisotropic model detailed in equations (13) and (15), and find that the differences are very similar in magnitude, as is expected by the small difference () between the momentumtransfer crosssection that enters the inelastic collision operator (equation (II.2.2)) for isotropic scattering, and the anisotropic form.
For each of the cases tested, our calculations demonstrate that significant differences can appear when various approximations are made, or detailed balance is neglected. Of particular interest in this study is anisotropic scattering in the lowthreshold inelastic channel for our model system, where differences up to 11% in the transport coefficients were calculated from the inelasticallyisotropic model.
iii.2 Electron transport in molecular nitrogen
iii.2.1 Crosssection set
The set of crosssections utilised throughout this work were obtained from the v8.97 Magboltz database tabulated on LXCat (Biagi, 2012). Other databases of scattering crosssections for electron swarms in molecular nitrogen are available (e.g., see (Pitchford et al., 2017)), but our focus here is to illustrate the damping effects of temperature on NDC in a real gas rather than an analysis of the selfconsistency of the different sets. The crosssection set (Biagi, 2012) details elastic collisions, 15 individual vibrational states, 29 electronic state collisions and a single ionisation crosssection, with rotational excitations to be included as in the Magboltz source code (Biagi, ), outlined below. We also note that the elastic momentumtransfer crosssection tabulated on LXCat is more sparse at lower energies than in the source work of Itikawa (Itikawa, 2006), so the tabulation from Itikawa at lower energies is preferred.
The LXCat Biagi database at present does not include rotational crosssections, and although other databases do include a description of rotational collisions, the rotations of the source work (Magboltz) are preferred for their consistency with the elastic momentumtransfer crosssection. The rotational excitations utilised in the Magboltz v8.97 source code (Biagi, ), to be included with the LXCat tabulation for low field simulations, were calculated using the Gerjuoy and Stein treatment of the atoms as point quadrupoles with the Born approximation (Gerjuoy and Stein, 1955a, b). Rotational crosssections for the transitions were considered, using the values for the quadrupole moment constant in units of , where is the Bohr radius, and the rotational constant eV (Itikawa, 2006; Biagi, ). In Magboltz, an enhancement of the crosssection magnitude in the resonance region between 1.2 eV and 5.3 eV was explicitly included in the rotational crosssections. This is also included here, along with scaling of the rotational crosssections above 5 eV to fall at the same rate as the elastic momentumtransfer crosssection. For the temperatures considered in this work, sufficient convergence in the calculated transport coefficients is obtained using 40 rotational transition crosssections, and we adopt no fewer here. Note that the Born approximation, Gerjuoy and Stein, rotational crosssections are arguably valid up to 0.6 eV (Gerjuoy and Stein, 1955a). For the transport coefficients considered in this work, up to 400 Td, the rotational crosssections are each calculated up to 400 eV. At greater than approximately 20 Td, however, when vibrational processes are active, the power (figure 10) and momentumtransfer rates (not shown) indicate that rotational collisions, both inelastic and superelastic terms, are considerably less important than for the other processes.
We also consider the energy sharing fraction between the two postcollision electrons resulting from the ionising collisions, taken here to be equally shared between the scattered and ejected electrons. At 360 Td, the highest reduced electric field measured in our experiment (see figure 9 and Appendix B), we find a difference of less than 0.4% and 0.6% between the drift velocities and mean energies, respectively, from the 50%50% sharing fraction and 1%99% sharing fraction. Comparing the 50%50% sharing fraction results with those for allfractions being equiprobable, a less than 0.02% difference is calculated between the drift velocities, and less than 0.3% between the respective mean energies. The size of these differences is not unexpected given that the power transfer from the ionisation channel is of a similar magnitude to the (individual) electronic state excitations at the highest reduced electric field considered, as shown in figures 10 and 11.
iii.2.2 Temperature dependence of the electron transport properties in N — measured and calculated
In figure 9 and Appendix B we present our experimental values of the drift velocity and the Townsend ionisation coefficient for electrons in molecular nitrogen, as a function of the reduced electric field in units of the Townsend (1 Td = ). The experimental drift velocities measured at 293 K and over the range 0.65–360 Td are in reasonable agreement, but tend to lie below, the other available experimental data (by at most 10%, with the exception of the Wedding et al. (Wedding et al., 1985), Roznerski (Roznerski, 1996; Snelson and Lucas, 1975), and Kelly (Kelly, 1990; Campbell et al., 2001) data at the higher where that difference increases up to 19%). Our measurements also lie below our calculations at 293 K, with differences of less than 2.8% over the range of reduced electric fields measured, which is somewhat larger the overall uncertainty of .
Our experimental measurements of the Townsend ionisation coefficient compare reasonably well with the other available experimental measurements; within 5% of the measurements of DeBitetto and Fisher (DeBitetto and Fisher, 1956), within 14% of the HernándezÁvila et al. measurements (HernándezÁvila et al., 2004a, b) and those from Cookson et al. (Cookson et al., 1966) and Kelly (Kelly, 1990; Campbell et al., 2001), and generally underestimating the other experimental measurements shown in figure 9 by less than 50%, with the exception of the measurements of Haydon and Williams (Haydon and Williams, 1976), and some of the Wedding et al. (Wedding et al., 1985) measurements, although the bulk of these data lie within 5% of our present measurements. Comparison between our experimentally measured Townsend ionisation coefficient and our calculations showed similar discrepancies. Namely, our measurements are lower than our calculated values over the range of reduced electric fields 120–360 Td, with a difference around 55% at the lower fields where the coefficient is rapidly rising, and decreasing to a difference of around 14% at 360 Td. These differences are larger than the experimental overall error bars of . Compared with the other available measurements, our calculations tend to overestimate the experiment below 200 Td by generally less than 60%, decreasing with increasing field, to be within 15% at the higher reduced electric fields.
Comparing our transport coefficients calculated at 293 K with our calculations at 300 K, they show an up to 3% variation, decreasing below 0.1% above 1 Td, in the drift velocity and mean energies, and a less than 0.01% difference in the Townsend ionisation coefficient above 100 Td, increasing to 1.25% at 60 Td.
The observed discrepancies between our calculated and measured transport coefficients are addressed in the following section.
iii.2.3 Approximation effects: Twoterm approximation, drift velocity definition, and higher order superelastic processes
Before our discussion of the effect of temperature and superelastic populations on NDC in N, we digress to consider the effect of some of the approximations associated with calculating transport coefficients for electron swarms in N. This subsection details the effect on the calculated transport coefficients of the different definitions of the drift velocity, the twoterm approximation, and the treatment of superelastic collisions.
In consideration of the differences between our calculations and the experimental data above, we have investigated the effect of on the calculated drift velocity and the Townsend ionisation coefficient, specifically using a twoterm solution. The limitations of the twoterm approximation have been discussed in detail previously (e.g., (White et al., 2003)) and our calculations here illustrate some of these differences. Figure 12 shows our results using twoterm and multiterm Boltzmann calculations, where the agreement with our experimental measurements is improved by using a twoterm solution. The errors decrease from less than 2.8% and 14–56% difference, between our experimental drift velocity and the Townsend ionisation coefficient for our multiterm calculations, respectively, to less than 2.6% and 8–28% for the twoterm calculation results. While this may initially appear to be counterintuitive, in fact it simply reflects that the Biagi (Biagi, 2012, ) crosssection database we used was originally engineered for application with a twoterm code to reproduce a selection of the available measured transport coefficients. For all calculations other than in this figure, the results presented are from multiterm calculations.
We also note that the neglect of recoil in the inelastic channel in our solution may have a similar impact on our N calculations as for our model crosssection. The mass ratio and rotational threshold in the model are similar to those for N, and the power transfer rates show similar behaviour to the model calculations. This observation may be able to account for some of the underestimation of our calculated transport coefficients when compared with those from the present experiment.
Transport coefficients are dependent on how the experimental current trace is analysed (Robson, 1991; White et al., tion), so in figure 12 we compare the calculated flux, bulk, and steadystate Townsend drift velocities to the drift velocity extracted from our pulsed Townsend experiment. It is expected that the differences between the various possible drift velocities increase with increasing reduced electric field, with the particle nonconserving ionisation channel increasing in importance (as shown in the power transfer rates in figures 10 and 11). At the highest measured reduced electric field of our new experimental data, 360 Td, the difference between the flux and SST drift velocities is 2.5%, while a 10% difference is calculated between the bulk and flux drift velocities at this field.
For our calculated Townsend ionisation coefficient we find less than a 0.1% difference between the coefficient calculated from a second order approximation to the Townsend ionisation coefficient using the bulk timeofflight coefficients (given in equation 11), and the direct calculation of the coefficient using our steadystate Townsend simulation over the range of our experimental values, as shown in the lower pane of figure 12.
The neglect of deexcitation processes in the terms of the inelastic collision operator, but inclusion in the term, is considered in N to assess their effect on the drift velocity and the Townsend ionisation coefficient. Unlike with our model calculations in section III.1.3, with only one excitation channel, in N the scaling of the groundstate excitations, when modifying state populations, requires more consideration as the number of rotational channels significantly populated changes with the temperature of the neutrals. There are multiple scalings of the groundstate equation that may be considered (for example, with nonzero temperature or at 0 K, or with or without degeneracy considered), so we have assessed the two extreme possibilities combined with setting the superelastic population to zero for . We first consider using the proper groundstate density for , calculated using the neutral temperature, and secondly with no scaling at all (effectively ):

When the density of neutrals in the ground state are calculated according to the temperature of the neutrals, the differences between our standard 293 K calculations of the drift velocity and the Townsend ionisation coefficient change by less than 4% and 0.1–31%, respectively. The differences between both coefficients decreases with increasing reduced electric field, as the higher threshold processes with lower excitedstate densities (being neglected) starting to dominate.

For the extreme case of no temperaturedependent scaling or degeneracy included in the equation, when setting , this is not equivalent to the case considered in our model calculations with 0 K in the equation of the inelastic operator. In this case, all of the rotational processes would be weighted equally, so the resulting transport would be influenced by the number of rotational crosssections included in the set. In our current set of more than 40 individual rotational processes, the differences between our standard calculations and this modified set, for the drift velocity and the Townsend ionisation coefficient are up to 60% and 8–100%, respectively. The much higher number of neutrals in the ground state for each excitation channel results in higher momentum exchange with electrons with energies reduced to , resulting in a reduced mean energy and drift velocity, and a reduced Townsend ionisation coefficient as sampling of the ionisation crosssection is delayed to higher reduced electric fields.
For both of these extreme cases, the neglect of superelastic processes has an important impact on the calculated transport coefficients, particularly the Townsend ionisation coefficient, when considering the 0.1% accuracy required in swarm calculations, and even the 10% error acceptable in plasma applications (White et al., 2003).
iii.2.4 NDC in N
In molecular nitrogen, an NDC region is present in the measurements of both Pack and Phelps (Pack and Phelps, 1961a) and Lowke (Lowke, 1963) at 77 K and 77.6 K, respectively. Our calculations and those of Petrović et al. (Petrović et al., 1984) are in good agreement with the experimental measurements of these authors. Petrović et al. (Petrović et al., 1984) also presented drift velocity calculations at 293 K with and without superelastic collisions, that our calculations are in similarly good agreement with. An increase in the drift velocity and the appearance of an NDC region occurs in the absence of superelastic processes. Highlighted in their work is the importance of superelastic collisions for both the lowthreshold rotational and vibrational excitation channels. As with our model crosssection, the temperature dependence of the excitation state populations in N, and the reduced nett energy transfer due to inelastic collisions, compared to the 0 K calculations, is responsible for the absence of NDC at 195 K and 293 K. At 77 K, the decreased population of the excited states for rotational and vibrational excitations results in a higher nett energy transfer rate in those channels, and a corresponding decrease in the drift velocity with increasing reduced electric field. In figures 10 and 11 the power transfer rates for each collision type are given for each of the temperatures considered. The explicit power transfer due to the deexcitation processes (in particular rotational excitations, here grouped into a single line for the figures only) illustrate that the increased contribution of the superelastic processes at the higher temperatures is responsible for the decreased range and eventual disappearance of the NDC region, where the increased mean energy of the swarm changes the sampled region of the elastic collision frequency, resulting in an increased drift velocity. The temperature dependence of the NDC region has been observed previously through the vibrational channel temperatures in mixtures of 99% argon and 1% N (Dyatko et al., 2010), and our present discussion highlights the same dependence occurring in a pure gas.
Iv Concluding Remarks
In this study we have investigated the temperature dependence of NDC using a simple model system, alongside that of N. The power transfer rates in the elastic, inelastic and superelastic channels show the damping effect of the deexcitation processes on the range of NDC. With increasing temperatures, the higher proportion of neutral background gas particles in excited states increases the mean energy and subsequently suppresses the NDC region that arises from the increasing elastic momentumtransfer crosssection compared with the (decreasing importance of the) inelastic channels at those fields.
To assess the impact of superelastic collisional processes on NDC, we employed some model (although unphysical) cases, with temperature dependence during elastic, excitation and deexcitation processes manipulated, illustrating the importance of the deexcitation process to the transport coefficients at low reduced electric fields. These systems also isolated the physical processes responsible for NDC, with the energy gained by the electron swarm from the deexcitation channel reducing the range of or eliminating NDC altogether.
We have also presented calculations of Robson’s (Robson, 1984) criterion for the presence of NDC using the rate of change of the ratio of the nett energy exchange of inelastic to elastic collisions, derived using momentum transfer theory. That criterion predicts well the region of NDC in all of the model cases considered, as well as the temperaturedependent NDC region present in N, using only a knowledge of the collision frequencies.
The effect of anisotropic scattering, for verylow threshold inelastic processes on the transport coefficients, was assessed using a model crosssection to replicate the forwardpeaked nature of rotational collisions. The inclusion of an inelastic momentumtransfer crosssection results in a 5–11% increase in the drift velocity, and between a 4–9% change in the mean energy of the electron swarm for the temperatures considered in this work.
In the FrostPhelps differential finite difference form of the inelastic collision operator utilised in this work, the representation of inelastic collisions is truncated at zeroth order in the mass ratio, neglecting the recoil of the neutral particle during an inelastic collision. The effect of this assumption had been assessed previously and found to have less than a 0.1% impact on the calculated transport coefficients for electrons in H (White et al., 2002). For the model crosssection considered in this work, however, the inelastic threshold is more than 20 times lower than the lowest rotational threshold in H, and the impact of the truncation of the mass ratio for inelastics was found to have a greater influence on the transport coefficients. At 0 K recoil accounts for a less than 2% change in the drift velocity and mean energy, but this difference increased to over 6% at room temperature. To derive the next terms in the mass ratio expansion for the FrostPhelps inelastic operator was beyond the scope of the present work, however should be considered when adjusting crosssections derived from swarm transport measurements for processes with very low thresholds (for example, the derived vibrational crosssections for H (White et al., 2002, 2007)).
Finally, we have reported experimental measurements of the drift velocity and the Townsend ionisation coefficient for electron swarms in gaseous N using a pulsed Townsend apparatus. Comparison of our measurements with some of the other available experimental measurements shows reasonable agreement, generally within 10% for the drift velocity, and 5–50% for the Townsend ionisation coefficient. We have also presented our calculations using a multiterm Boltzmann solution, that well reproduce experimental drift velocities and the Townsend ionisation coefficients at 293 K, but tend to somewhat overestimate our current experimental measurements. At lower temperatures, we also tend to underestimate the experimental drift velocities at the lower reduced electric fields, however addressing this by modification of the N crosssections we employed through a swarm analysis was outside the scope of this work. Rather, our calculations were used to illustrate the physical processes associated with NDC and the effect of temperature on its appearance or absence, with the same dependence on superelastic populations found as in our model calculations.
Acknowledgements.
The authors would like to thank the Australian Research Council through its Discovery Program (DP180101655) for financial support. The experimental study was supported by UNAMPAPIIT IN 108417 and Conacyt project 240073. The technical assistance of A. Bustos, G. Bustos and H. Figueroa is greatly acknowledged. SD would like to acknowledge MPNTR projects ON171037 and III41011 for support.Appendix A Numerical Considerations
The theory and solution techniques employed in this study have been systematically benchmarked against independent Monte Carlo and kinetic theory solutions. Using model systems, each of the collisional processes have been validated by comparing against existing solutions for the hard sphere and Maxwell’s models, Reid’s ramp and anisotropy models (Reid, 1979), and attachment and ionisation models (Ness and Robson, 1986). When we investigate these model systems outside of the reduced electric field and temperature ranges generally considered, we sometimes find problematic behaviour in their solutions and resulting transport coefficients. For some real gases we are also able to provoke similar behaviour, where the electron energy distribution function contains a contribution from a solution that is not part of the physical solution we expect to extract, as illustrated below. Generally these problematic regions are outside standard swarm experimental configurations, but warrant further investigation nonetheless.
To assess the effect of the nonphysical contribution on our calculations, for all of our results here we compare with an independent Monte Carlo solution, discussed in section II.5, and find that for the conditions required for calculations involving real gases, and the model system employed here, we reproduce the Monte Carlo results, with the exception of the NDC region of our model system where recoil of the neutral during inelastic collisions becomes important, as discussed in section III.1. As part of our investigation into the origin of these problematic solutions, we have trialed different boundary conditions from those of Winkler and collaborators (Loffhagen and Winkler, 1996) that we usually employ. The generalised eigenvalue method utilised for our timeofflight and steadystate Townsend solutions can be solved using inbuilt Matlab functions, or a benchmarked inverse power method, where we enforce strict convergence criteria. Using a high minimum number of iterations in the inverse power method solver, improvements in the distribution function are obtained, but the contributions from the nonphysical solution remain.
We illustrate using N that a discrepancy exists in our simulations below the experimental values. For timeofflight calculations, when solved as a generalised eigenvalue problem, the lowest temporal eigenvalue is equivalent to the nett rate coefficient, but is not explicitly used in calculations of the transport coefficients, given in section II.4. For the steadystate Townsend simulations, however, the lowest spatial eigenvalue represents the Townsend ionisation coefficient extracted from the experimental measurements. For any configuration with conservative collisions only, both spatial and temporal eigenvalues should be zero, and we use this, along with the equivalence of the nett rate coefficient and temporal eigenvalue, as a consistency check for our solution.
Depicted in figure 14 for N is the Townsend ionisation coefficient calculated using these equivalent methods. In the low field region (below our present experimental measurements) the ionisation crosssection is sampled by the tail region of the distribution function only, so we would expect to see the Townsend ionisation coefficient decreasing with decreasing reduced electric field. Using the bulk transport coefficients from a timeofflight simulation in a second order approximation to the spatial rate coefficient (given in equation 11), our calculations well reproduce the experimental results, labelled ‘TOF , standard’ in figure 14.
With the temporal eigenvalue used in place of the equivalent nett rate coefficient ( in equation 11), labelled ‘TOF e’val, standard’ in figure 14, in the low field region our calculations produce nonzero values for the Townsend ionisation coefficient, showing an inconsistency in our timeofflight simulation.
In our steadystate Townsend simulation (labelled ‘SST, standard’ in figure 14), where the lowest eigenvalue corresponds directly to the spatial ionisation rate coefficient, we observe the exact same nonzero contributions at reduced electric fields below where ionising collisions should be contributing. The consistency of this nonzero contribution between the two different simulations suggests a leak of number density stemming from the same numerical source.
To address these issues, in the following subsections we consider the numerical representation of the collision operators, where using an alternative description of the collision operators results in electron energy distribution functions where the noise in the tail region is suppressed to satisfactory higher energies/lower magnitudes. In regions where the distribution function does not drop a satisfactory amount, all of our results are calculated using these alternative collision operators and are compared with an equivalent Monte Carlo solution, see section II.5. Under these circumstances we find very good agreement between the distributions and resulting transport coefficients. This is illustrated in N by the Townsend ionisation coefficient calculations labelled ‘Conservative’ in figure 14, where using the timeofflight coefficients with either or the equivalent temporal eigenvalue, and our steadystate Townsend code reproduce the expected values. We note that each of the modifications proposed below reproduce the standard model benchmark systems.
a.1 Representation of the derivative in the elastic collision operator
For elastic collisions with nonstationary neutrals, the collision operator has a single and double derivative term in the expression. Using the collision operator when the first derivative term is expanded with the product rule meets all of the existing benchmarks. Similarly, if we now numerically represent the collision operator as , again all benchmarks are met, and for some situations the contribution of the noise in the tail region of the electron energy distribution is delayed to higher energies, allowing the solution to capture a sufficient energy range of the electron swarm. This is illustrated by an example using Maxwell’s elastic model crosssection at 293 K with Td, with the resulting and terms of the distribution function shown in figure 15.
This representation is equivalent to a finite volume approach, where a flux of electron density is moved between energy grid elements, with conservation of this density enforced. For our illustration in N in figure 14, using the timeofflight coefficients and temporal eigenvalue in place of , this conservation of density reduces the Townsend ionisation coefficient calculated at reduced electric fields below where the ionisation channel has a significant impact on transport, however a nonzero contribution remains, to be addressed in the following section that considers the representation of inelastic collisions.
Similar considerations for the representation of the field term may be taken, representing the derivatives together as other authors employ (e.g., (Loffhagen et al., 2002)), or expanding using the product rule (as in reference (Trunec et al., 2006)), but we have found very little difference in the resulting distributions or transport coefficients.
a.2 Conservation of number density in inelastic collisions
To address the inconsistency between the eigenvalue and nett rate coefficient, highlighted for N in figure 14 above, we reformulate the inelastic collision operators to force conservation of electron density in a finitevolume style. To calculate the change in the distribution function , due to an inelastic collisional process , an expression is needed for the change of the number density over time , in terms of known quantities, , and calculable quantities, or .
For the flux of particles in and out of each energy bin element to be conserved, to move (or add for ionising collisions) particles between pre and postcollision positions in the energy grid, a movement matrix is applied to move the density of particles in each volume element surrounding the solution energy grid . The general form for the change of the distribution function due to an inelastic collision is then given by , where is the neutral number density, scales the crosssection in m to , and are the bin widths (i.e. ). The energy dependent quantities to the left scale the conserved number density (contained in ) per energy element per unit time, to the distribution function modified by the collision . Here, the conserved quantity is the number density scaled by the collision frequency for the process , taken to be .
The general form of the movement matrix calculates the density of particles from the prescattered energy element , with associated bin , that move into each of the energy bins in the postcollision scattering region , and varies for inelastic, superelastic and ionising collisions. For an expression for the flux of the density of particles in and out of each bin, we consider the general case where the boundaries of the pre or postscattered energy bin lie within (not at the edges of) a volume element surrounding a point in the solution grid (i.e. for an element of the solution grid to be scattered, and postscattered energy bin with left and right boundaries , taken to lie within a grid element , and ). The number density between a left and right boundary can be calculated using the density at the bin edges, and here we choose to linearly interpolate the distribution function. For a bin , the density of particles, scaled by the collision frequency for the process , taken to be , is given by . When is assumed to be linearly spread across each bin width, the expression for at some energy is given by: , so that