Influence of Elastic Scattering on Electron Swarm Distribution in Electrified Gases.

The propagation of energetic electrons through air is one key component in the generation of high-energy atmospheric phenomena such as lightning-generated X-ray bursts, terrestrial gamma ray flashes (TGFs), and gamma ray glows. We show here that models for this propagation can be considerably affected by the parameterization of the differential cross section of elastic scattering of electrons on the molecular components of air. We assess existing parameterizations and propose a more accurate one that builds upon the most up-to-date measurements. Then we conclude that by overweighting the forward scattering probability, previous works may have overestimated the production of runaway electrons under high electric fields close to the thermal runaway threshold. ©2020. The Authors.

Introduction

The discovery that in a thunderstorm environment electrons are sometimes accelerated to relativistic energies (Dwyer et al., 2012) has prompted a renewed interest in the propagation of energetic particles through air. There is consensus that electron thermal runaway is the key process behind the X‐ray emissions from lightning (Dwyer et al., 2011; Moore et al., 2001) and in long sparks in the laboratory (Dwyer et al., 2005; Kochkin et al., 2012, 2015, 2016; Montanyà et al., 2015; Rahman et al., 2008; da Silva et al., 2017) and possibly it also plays a role in the production of terrestrial gamma‐ray flashes (Fishman et al., 1994). Thermal runaway is possible because below a few megaelectronvolts the rate of energy loss by an electron in air peaks around 300 keV/cm for an electron energy close to 200 eV (Dwyer et al., 2012; Gurevich et al., 1999). This implies that under an electric field above 300 kV/cm on average all electrons accelerate, even those in the bulk of the energy distribution function, and thus their maximum energy is only limited by the total available electrostatic potential or by radiative deceleration (bremsstrahlung) if they reach a few megaelectronvolts (Köhn et al., 2017; Xu et al., 2012).

There are several competing hypotheses about the mechanism that produces fields above 300 kV/cm and initiates thermal runaway in a discharge. Many authors have investigated streamer head‐on collisions (Babich & Bochkov, 2017; Cooray et al., 2009; Ihaddadene & Celestin, 2015; Köhn et al., 2017; Luque, 2017; Lehtinen & stgaard, 2018), whereas others have proposed the propagation of ionization waves along preexisting channels (Babich et al., 2015) or a reduction in air density, which allows for thermal runaway under proportionally lower electric fields (Köhn et al., 2018). As none of these models has gained a consensus among the research community, it is presently not clear whether any of them can explain the existing observations. Therefore, we need either a better macroscopic modeling of the streamers within a corona or a better microscopic modeling of the dynamics of single electrons in air.

In this work we consider the second possibility and focus on the statistical distribution of scattering angles after an electron interacts elastically with a molecule. We discuss existing measurements and approximations to this distribution, which is the angular part of the differential scattering cross section.

The paper is divided as follows: In section 2 we review the literature about the differential cross section (DCS) of elastic electron scattering in air and build a new, more accurate parameterization suitable for numerical models. Section 3 describes a numerical Monte Carlo code for the propagation of electrons in air. We use this code to show, in section 4, the effect of different parameterizations of the DCS on the energy distribution function of electrons under a high external electric field, underlining the significant effect in the tail of the distribution. Finally, we offer some conclusions and recommendations in section 5.

Elastic DCS Representation

Figure 1 sketches the geometry of an electron scattered off a molecule. The outgoing direction of the electron is defined by the polar angle of deviation θ and the azimuthal angle φ. If the electron advances with a speed v (kinetic energy ε) within a uniform number density n of molecules, the (differential) probability per unit time that the electron scatters into a differentially small solid angle , is nv(dσ/dΩ)dΩ, where dσ/dΩ is called single‐scattering DCS. Here we are concerned with elastic scattering, where the internal energy of the molecule is not altered by the collision. In that case, due to conservation of momentum and energy, the electron's initial kinetic energy ϵ is partly transferred to the target molecule so that its final energy ϵ is determined uniquely (Canto & Hussein, 2013b) by the scattering angle θ and the mass ratio μ of the electron to the molecule:

Figure 1

Sketch of electron elastic scattering by a target: The outcome of the collision is described through the polar deviation angle θ with respect to the incident direction, and an azimuthal angle φ that obeys a uniform distribution corresponding to the cylindrical symmetry assumed. The distribution in θ is described by the differential cross section dσ/dΩ.

Across the range of energies covered, from a few electronvolts to kiloelectronvolts, the qualitative shape of the DCS varies considerably. We will often refer to three ranges of electron energy:

low‐energy: within a few sigmas from the average kinetic energy (up to 20 eV and below),

middle‐energy (or intermediary): where backscattering is still observed in the DCS (from 30 eV up to 200 eV),

high‐energy: passed the friction curve peak where electrons have high chances of becoming runaways (from 500 eV).

Note that this classification is not exhaustive and that there are intermediate regimes between these ranges.

Figure 2 shows the experimentally measured DCS for electrons scattering on molecular nitrogen for three energies. At 1.92 eV the scattering is mostly isotropic from where it becomes gradually replaced by forward scattering as the energy increases. From 50 eV one also notices an exponential decrease of the DCS at very small angles, attributable to cloud polarization (Daimon et al., 1982; Herrmann et al., 1976). In the intermediary range, there is also a residual backscattering effect that vanishes with increasing energy. At 500 eV, we can observe already undulations at large angles interpreted as coherent interference between wave functions scattered about displaced individual atoms (Hayashi & Kuchitsu, 1976a).

Figure 2

DCS data for N2 (in m2 per solid angle) in three different energy regimes: low energies yield an isotropic‐like distribution (relative error ∼ 7%), middle energies present forward and back scattering peaks (relative error ∼8%), and high energies are marked by forward scattering (relative error ∼10%).

In our study we gathered the most recent experimental measurements available of elastic DCS for N2 and O2, some of which were already summarized by Trajmar et al. (1983) and Brunger and Buckman (2002). Although very useful, the data are restricted for practical reasons to energies below 1 keV and angles ranging from roughly 5° to 160° where relative errors range in most cases between 5% and 15%. Both extremes of low and high angle values usually have larger uncertainties than those at middle‐range angles. All the DCS were fitted by ordinary least squares where datapoints were weighted by the inverse of their squared errors. In some of the figures containing experimental DCS those error bars will not be displayed for matters of clarity.

To overcome those limitations and acquire more dense and extended DCS, one must resort to model calculations. The theory of electron‐molecule scattering is well developed (see, e.g., Burke et al., 1977; Hayashi & Kuchitsu, 1976b; Sun et al., 1995) and provides different methods to calculate the DCS in various energy ranges. Nevertheless, there are two major inconveniences against using those models in Monte Carlo codes. First, no single model can accurately and efficiently compute DCS over the energy range covered here (eV → keV). Different approximations would be used including the Born approximation beyond 1 keV, the independent atom approximation at several hundred eV (Hayashi & Kuchitsu, 1976b), and several exchange and correlation‐polarization potentials below 100 eV (Kohl et al., 1981). The methods used for solving those approximations are not the same. Second, those methods might be time consuming for one given electron energy. Hence, the tabulated data may tend to be sparser in energy than the desired resolution.

In practice, most often the DCS is fitted with a (semi)empirical formula. In some very simple cases, the formula can be inverted to sample analytically the scattering angle. Otherwise, the formula can be used to generate a 2‐D array of DCS values at the desired energy‐angle resolution. Starting from a few tens of electronvolts up to tens of MeV beyond which relativistic corrections may apply, most approximations build upon the so‐called Screened Rutherford (SR) cross section, which estimates the DCS for the scattering of an electron around an atomic nucleus: where m and v are the electron mass and velocity (kinetic momentum p=mv), the charge of the target is Ze, ε 0 is the vacuum permittivity, and η is a so‐called screening factor. This parameter accounts for the screening of the charge in the nucleus due to surrounding electrons. By reducing the strength of the interaction at large distances, the effect of η is to reduce the scattering at small angles θ. In the limit η→0, the expression (2) reduces to the diverging Rutherford (R) scattering. The emergence of the SR comes from applying the first Born approximation on a Yukawa‐type potential representing the atom. The latter potential was originally derived by Wentzel (1926); more details are given by, for example, Sakurai (1994).

The SR cross section is very appealing for its simplicity and transparency. For this reason, in the realm of electron transport in Monte Carlo codes (Berger et al., 1970; Chanrion & Neubert, 2010; Jackman & Green, 1979; Köhn & Ebert, 2014; Köhn et al., 2017; Moss et al., 2006), the DCS used came from various analytical models (Ivanov et al., 1977; Jackman & Green, 1979; Murphy, 1988; Moss et al., 2006; Porter & Jump, 1978), all of which were, at least in part, built upon the SR cross section; for either scattering about molecular nitrogen or oxygen. Following this trend, we have revisited those formulae and built our own DCS fit for Monte Carlo‐based simulations.

It is our objective to challenge whether the DCS can accurately be extrapolated with the SR for the application of electron transport in very high electric fields (a small fraction of the critical value for thermal runaway). To begin with, we will compare the shapes of various DCS functions and relate these with the corresponding probabilities of scattering an electron. Then, we will show to what extent macroscopic properties differ when each of the aforementioned analytical DCS are used in simulations. Since the earliest DCS versions were criticized by the latest we will only present and recur the following three. Although those models wind back to several decades ago, they are still in use as benchmarks for Monte Carlo codes of electron transport in recent years.

Jackman and Green

In a study about electron degradation in atmospheric gases, Jackman and Green (1979) presented and compared an improved DCS to the basic SR. They took into account the experimental evidence (Bromberg, 1969; Daimon et al., 1982; Herrmann et al., 1976) that the DCS for small angles at higher energies follows a negative exponential trend in the angle θ due to the polarization of the atomic electron cloud induced by the incident electron (Herrmann et al., 1976). Also, to include backscattering below 200 eV for N2 and O2, they added a term proportional to the SR cross section with the substitution θ→π−θ (i.e., ). This term adds a peak at high scattering angles; we will refer to such term as “reverse screened Rutherford” (RSR) from now on. The full expression reads

This formula is normalized to the absolute cross‐section σ (ε) and comprises several parameters that are functions of energy: with the numerical values listed in Table 1. Note that there is a sign error in their original formula for the normalization of the exponential term in (3). The authors mention the applicability of this formula only for N2; here we apply it to O2 as well.

Table 1

Jackman and Green Parameters for DCS in (3)

f 11	f 12	f 13	a 1	a 2	a 3
100 eV	0.84	1.92	1 eV	0.43	−0.29
f 21	f 22	f 23	c 1	c 2	c 3
10 eV	0.51	0.87	12 eV	1.27	0.27

Their formula (3) is reported to be valid and verified only down to 30 eV. It breaks at 12 eV because of the pole introduced when η 2(ε) reaches a negative value. Jackman and Green (1979) used a multiple scattering formula for lower energies, despite most multiple scattering theories being based on small angle scattering (Goudsmit & Saunderson, 1940b; Moliere, 1948; Nigam et al., 1959; Scott, 1963) and thus more appropriate for higher energies. In our tests, we chose to use their formula down to 15 eV and cut it with an isotropic distribution below that threshold. This choice was made for convenience; as their fit at 15 eV is still a better approximation than the isotropic distribution.

Another unusual feature of (3) is that the screened parameter η 1 is fixed to a constant value; small‐angle scattering is controlled by an energy‐dependent balance between the first (exponential) term and the second (SR) term. As this formula is not founded on a physically consistent extrapolation of the atomic screening, it is restricted to energies below a few kiloelectronvolts. In the simulations performed by Jackman and Green (1979) electron energies did not exceed 5 keV.

Murphy

Similarly to Jackman and Green (1979), Murphy (1988) uses an SR for forward scattering and a RSR for backscattering: His normalized DCS is with

In this expression the factor N(ε) is interpolated linearly from precomputed values (Murphy, 1988). As in the previous DCS case, the formulae (9), (10), (11) above originally apply only to N2; here we use them also for O2.

As opposed to (3), (9) converges naturally toward an isotropic distribution at low energies due to the unbounded values taken by the screening parameters η 1,2. The main shortcoming of (9) is that it is based on a single experimental set of measurements by Shyn et al. (1972), which is limited to electron energies from 5 to 90 eV.

Moss et al.

Moss et al. (2006) used a more extended database of DCS, on which the present work is partly based. For energies below 500 eV for N2 and 1 keV for O2, they directly computed the DCS from the data; supposedly through linear interpolation. For higher energies, they used a modified SR (MSR (12)) (numbered (16) in their paper) with the screening parameter off bracket in the denominator. This modification steepens the demarcation line between forward and large‐angle scattering: where η=2ε 1/ε and ε 1=4eV.

Present Fit

In all the formulae presented above the SR DCS dominates at high energy. However, each of these formulae has its weaknesses. For instance, (3) cuts sharply the backscattering at 200 eV and presents a large discontinuity at the threshold of isotropic scattering, whatever its value; (9) is subject to first derivative discontinuities induced by the linear interpolation of N(ε); and (12) was fitted and snapped to N2 at 500 eV without any transition. Finally, all the formulae above apply strictly to N2: their application to O2 is implicit.

Aware of the varying shape of the DCS with energy (cf. Figure 2), we expanded the DCS into a traditional screened Rutherford (SR) term to account for the forward scattering, combined with a truncated series of Legendre polynomials (LP) in up to order L=3 to capture the remaining features due to a more complicated interaction with the scatterers and resulting mostly in large angle or backscattering. The main intuitive motivation behind using LP is that they correspond to axially symmetric spherical harmonics emerging from the assumed spherical symmetry of the average molecular scattering potential. For more information about the utilization of LP in DCS, the reader can refer to Goudsmit and Saunderson (1940a), Mott and Massey (1965), and Sakurai (1994). Our expression reads as follows:

For a given energy, this expression depends on two parameters for the SR term (the magnitude mag and the screening η) and on L+1 Legendre coefficients A . The magnitude parameter mag can be seen as correction applied to the coefficient in (2): , where m is of order unit. These parameters are adjusted to fit the database of experimental DCS values (evoked earlier and taken from  Trajmar et al., 1983; Brunger & Buckman, 2002) that cover the kinetic energy range from an ≈0.5 eV up to 1 keV. This range bridges the gap between a roughly isotropic scattering law for sub‐eV and a forward scattering law described by energetic electrons propagating in a screened Coulomb potential. Above 1 keV, the parameters are simply extrapolated and assumed applicable there too. The DCS is ensured to stay positive and follow the forward scattering trend exposed from roughly 200 eV to 1 keV.

The fitting proceeds in two major steps:

For each energy, the cross sections are fitted for all scattering angles (0° to 180° ), obtaining values for the parameters in (14).

For all energies, we then find expressions for these parameters as functions of the electron energy.

For energies below 15 eV, the DCS does not present a significant forward scattering peak and the SR term does not apply. In this case, the LP coefficients obtained at each energy were interpolated linearly and extended below the lower boundary of 0.5 eV by a fully isotropic scattering.

Beyond 15 eV, the full analytical expression devised for the DCS is the following: with and σ el(ε) is the absolute value of the cross section that is taken from the more precise database as given in Phelps and Pitchford (1985). For each element, the coefficients in (19) are displayed in Tables 2, 3, 4.

Table 2

Values of η Parameters

	A	B	C	D	K
N2	−1.56 ± −0.11	1.68 ± 0.16	2. ± 0.9	1.35 ± 0.09	−9.7 ± −0.43
O2	−1.07 ± −0.05	0.6 ± 0.53	1.4 ± 0.96	0.8 ± 0.22	−5.3 ± −0.36

Table 3

Values of r and for N 2

Al*	c l	m l	b l	d l	p l
1	1.6 ± 0.05	−6.8±1.3	1 ± 3	590 ± 370	−12.5 ± 8.6
2	2.6 ± 0.14	−20 ± 4.6	6.5 ± 1.4	2,200 ± 440	−41 ± 7.7
3	3.9 ± 0.36	−46. ± 14	−4.6 ± 3.2	5,900 ± 1340	−77 ± 15
r S	q 0=0.45±0.01		q 1=−5.10−4± 0.002

Table 4

Values of r and for O2

Al*	c l	m l	b l	d l	p l
1	1.38 ± 0.24	−1.3 ± 10	−110 ± 38	−15,000 ± 14,000	−120 ± 76
2	2.62 ± 0.14	−26 ± 5	14.7 ± 2.4	3,350 ± 480	−55.7 ± 7.4
3	3.94 ± 0.98	−70 ± 47	‐15 ± 10	18,600 ± 1,400	−138 ± 43
r S	q 0=0.49±0.04		q 1=−0.02±0.013

The newly introduced coefficients comprise the projection Sr (η) of the SR term in (14) onto the Legendre polynomial basis:

The details of the decomposition of SR into LP is given in Appendix A. We hereby introduce the following notation for the contribution of the SR to the l Legendre polynomial:

Since, just like in Murphy, 's 1988 case, we fit only the shape of the DCS and not its absolute value, the DCS is normalized so that all the reduced Legendre coefficients are defined as follows:

Figure 3 shows a capture of the functions fitting the Legendre coefficients for O2. With the exception of the zeroth order, we can observe for each graph in red solid lines what the contribution of the SR part (Sr (η)) to the total is. The first graph fits directly the ratio of the SR contribution to the norm ( , from equation (22)). A critical step in reassembling the analytic formulae for the Legendre coefficients and the SR is their relative magnitude. When tracing the values of (Figure 3, top left), the tendency at higher energies is not to vanish, but to the contrary: to maintain its magnitude erratically or even present a mild rise (Figure 3, top left). The choice for fitting is therefore difficult and we resigned to a simple logarithmic line (i.e., a power law). The standard deviation of these values relative to their least square linear regression is around 50%; which is very ominous. As we will see later (in Figure 8), this indicates that the SR cannot properly account for the scattering at large angles.

Figure 3

Normalized Legendre coefficients of the O2 DCS along energy. The first coefficient corresponds to r (23), the remaining coefficients represent the magnitude of extended Legendre coefficients (20) compared to the norm of the DCS. Each of the dots represents the coefficient of a Legendre polynomial that was fitted to experimental data at a certain energy. The dash‐dotted green curve represents fits (19) devised for points above 15 eV. In solid red: the SR decomposition into Legendre coefficients (see Appendix A).

Figure 8

Difference at high energies between experimental data and screened Rutherford. The Rutherford expression (R) is plotted along with the screened Rutherford (SR: equation (2)) and the modified screened Rutherford (MSR: equation (12)) in comparison with the DCS of N2 at 1 keV. The angular scale is given in for clarity.

Finally, because of the fact that does not vanish as the energy grows, the Legendre coefficients in (16) were cut off with a 1/(1+(ε/ε 0)) envelope so as to prevent spurious oscillations in the DCS at high energies.

An overall assessment of our analytical fit was conducted in Appendix B. There, we looked for two aspects: the (probabilistic) accuracy when used as a cumulative integral and the robustness/relevance of the parameterization. The robustness was interpreted as the standard deviation applying to the parameterization when the database was slightly deteriorated (reduction of datapoints). As can be deduced from Tables 3 and 4, some parameters are very ill defined due to the high standard deviation they have, and thus a better set of parameters should be sought in the future. The accuracy was measured by computing the average error (difference) between the exact and fitted cumulative integrals of the DCS at each energy (e.g., Figure 4). This error represents a probability shift in percentage for a population of electrons being scattered more (negative error) or less (positive error) than they should when they rely upon the linear interpolation of the database.

Figure 4

(left) DCS multiplied by for N2 at 60 eV, fitting by the SR and by our fit (16). (right) Cumulative integral of the DCS normalized. This illustrates the probability shift arising between an analytical fit and datapoints due to the extrapolation at small and larger angles. The difference between two curves relates to more (negative) or less (positive) probability of scattering of the first over the second curve. Here, we scatter less (above datapoints) than a linear interpolation of the data.

This probability shift is responsible for generating differences in the electron energy distribution under equilibrium in electrified air. We will see later how using only the SR as an approximation for the DCS at intermediate energies leads to very significant differences. We concluded that although the fit is very accurate and outruns any other of the analytical formulae covered, the parameterization has much room for improvement.

Ideally, the analytical formula in (16) can be used from 15 eV up to several hundreds of keV and can then be merged with the DCS relativistic formulation of Dalitz and Peierls (1951) for a Wentzel (Yukawa) potential. Eventually, a nuclear form factor can be introduced for very high energies (GeV) if one can find reliable models of nitrogen and oxygen nuclear densities (Olsen et al., 1964). It is nonetheless important to insist that beyond 1 keV, our extrapolation using the SR expression, equation (2), could not account for a proper reproduction of the DCS for the intermediate angles ( 40°) seen on experimental results. This is due to the nonvanishing LP contribution that could not be safely and reliably extrapolated beyond 1 keV.

Comparing DCS

All the presented formulae so far strive to provide a reliable representation of the DCS over a wide scope of electron energies. If we want to perform a pertinent comparison between those formulae, it is imperative to consider that the probability of deflecting in a polar angle θ must be weighed by due to the spherical symmetry of the problem sketched in Figure 1. This implies, as can be deduced from Figure 5, that extrapolations at very low or very large angles will only have a limited impact over the stochastic process of scattering compared to middle‐range angles close to 90°.

Figure 5

Extrapolation of N2 DCS at 4 eV: A linear interpolation (Linterp) is compared to various Legendre polynomial fits of different orders; the dashed lines represent the DCS fits weighed by .

For a reasonable range of extrapolated functions, the overall cumulative integral of the DCS as defined in section 3.1, equation (27), does not differ considerably as we can observe in Figure 5‐right. The crucial requirement is that the DCS does not display any aberrant behavior: such as to plunge into negative values in the case shown in Figure 5 of a LP extrapolation of order 8. When fitting a DCS, we strongly recommend to weight the data points with (or equivalently directly fit the data multiplied by ) as we did in this work.

With this first comment in mind, we compare (in Figure 6) the formulae presented so far in each of the three energy ranges presented in Figure 2. The graphs on the left show the DCS whereas those on the right show its cumulative integral (defined in (27)) that comprises the weighing by the factor. What is mainly of interest for the stochastic implementation of scattering is of course the cumulative integral, and so we can observe that various fitting techniques result in a different (shifted) probability distribution. Of interest, we also plotted what would the fits look like if we used Moss' formula (12) at 50 eV to be compared with the ordinary SR.

Figure 6

Comparison of several fits for N2 and O2 DCS and their cumulative integrals at various energies.

First, in the low‐energy part (Figure 6, top), Murphy's fit gets flatter to reproduce the isotropicity of low‐energy scattering. In the end, if we compare the efforts provided by a Legendre fit to a flat function, the difference is tenuous and probably even compensates at different (low) energies. As we move to the middle range (Figure 6, middle), Murphy's fit and this paper's fit correspond quite well to experimental data, Murphy did indeed use Shyn et al. (1972) data at 50 eV for N2 in his study. Jackman and Green's formula is slightly underestimating forward scattering because it strove to accommodate a broader range of energies into the same expression. As for Moss' formula, its essential difference comes from the fact that the low‐angle region is characterized by a plateau rather than a peak like in the case of the SR. Finally, at 500 eV (Figure 6, bottom), the differences build up, Murphy's fit, operated only from 5 to 90 eV is utterly out of shift and has a screening factor overestimation due to the scarcity of the data used (more detail is given in the next section). As for Moss's formula, it underestimates the scattering rate of midrange and high angles. Generally speaking, all formulae tend to underestimate the scattering rate at higher angles, and thus present a reduction bias of the backscattering. In our case, there appears a growing difficulty of fitting the backscattering remainder with a LP. Graphically, the LP part is supposed to fill in the gap between the blue (SR) curve and the data points. We will come back to this observation later on.

Screening Parameter

To better understand the implication of using the SR as the dominant forward scattering shape at high energies, we present a comparison of the screening parameters used by different authors. A logarithmic plot of η(ε) is given for Nitrogen and Oxygen (Figure 7). Values derived from a least squares fitting of the SR on experimental data are represented by dots. Then, various curves that try to represent η as a function of energy ε are plotted. We compare the previous results of Murphy (1988), Jackman and Green (1979), a simple linear regression (dotted line) and a hyperbolic fit (both in logarithmic scale), the screening parameter for multiple scattering of Nigam et al. (1959), and the theoretical formula for η corresponding to the approximation of the Thomas‐Fermi atomic potential by a screened Coulomb potential.

Figure 7

Logarithmic plots of the screening parameter η in (2) according to electron energy for N2 (top) and O2 (bottom). The additional experimental data used here come from Trajmar et al. (1983) and Brunger and Buckman (2002).

At first sight, the screening parameters seem to be adequately fitted by a simple power law (“Log‐Line”) as was performed by Murphy (1988). It can be demonstrated (Canto & Hussein, 2013a) that under the first Born approximation, there is an inverse square dependence of η∝ε −1 under nonrelativistic considerations. However, a closer inspection shows that there is a slight curvature trend to the log‐log plot (Figure 7). This is even more conspicuous for O2, where the low‐energy and high‐energy asymptotes remind of a hyperbola. There are at least three comments worth mentioning at this stage.

First, at lower energies, the influence of molecular structure prohibits the application of the first Born approximation from which the SR formula (2) is derived. Indeed, the eventual underlying screened potential is also perturbed with the more complex molecular potential, polarisation, exchange, etc. (Shimamura & Takayanagi, 1984). There is thus no reason for the screening parameter to correspond to its supposed theoretical value because the formula (first Born approximate) does not apply. For instance, the screening value might significantly differ from its current values in Figure 7 (top) if an exponential term is added to the fitting DCS, equation (14), and the SR is fitted on the corrected (subtracted) part. Thus, the values of η at intermediate energies, where polarization (and possibly exchange) effects are much stronger, should not be allowed to affect the trend displayed by those at much higher energies.

Second, most approximations for the exact atomic potential use a sum of Yukawa potentials (Ivanov et al., 1977; Moliere, 1947), instead of a single one. This means that the average screening parameter that encapsulates actually a squared sum of SR, can display a more complex dependence on energy than a simple inverse square law. This means that we are—voluntarily—using a generic formula (namely, the SR) onto the experimental database of N2 electron scattering, to observe the feasibility of encompassing the DCS; albeit it might actually comprise more terms.

Third, it has been shown (Bethe, 1953; Moliere, 1948; Nigam et al., 1959) that in the case of multiple scattering (with different targets) of electrons at high energies, the overall angular distribution f(θ,t) was found to depend on only two parameters: the thickness t of the path traversed, and a so‐called screening angle Ξ. The requirements for such a result were that the single‐scattering DCS respects or is close to the classical Rutherford formula at large angles (θ≫Ξ). The appellation of Ξ is due to the fact that (Bethe, 1953; Mott & Massey, 1965) when the single‐scattering law corresponds to SR (2). This essentially means that any complicated dependence of the DCS at small angles can be averaged through multiple scattering by using Ξ into a simple SR expression. This, however, is additionally conditioned by the requirement of using the SR as an absolute DCS and not a relative one. In other words, the use of Ξ in SR for Monte‐Carlo simulations implies that it is generally not correct to separate the DCS angular dependence and use another set of total elastic cross sections from a distinct database. This is because Ξ already encompasses and averages over the collision rate of a given (absolute) single‐scattering law. Consequently, the dependence of η(ε) upon energy can have several corrections to that of an inverse relation due to the fact that it may correspond to an averaging of a multiple scattered DCS. For comparison, we added to the single scattering η plot (Figure 7), the multiple scattering of Nigam et al. (1959) based on the second Born approximation of scattering about a screened Coulomb potential (Dalitz & Peierls, 1951).

Being aware that ideally, the power law of η(ε) should be very close to an inverse law (ε −1) for high energies—which might not apply to lower energies (as seen on Figure 7)—we used a geometric hyperbola in the logarithmic scale for fitting over the whole energy spectrum: which yields whereas Murphy's fit (10) would correspond to D=C=K=0 which yields a simple power law. The justification for this artifice comes from the fact that, upon realization that the SR is not a very good approximation of the DCS at intermediary energies, we wanted to separate the asymptotical behaviors at high energies from intermediate energies. One mathematical function known for comprising two asymptotes is the hyperbola. This hyperbola is applied on the logarithmic scale because the dependence of η on energy is a power law as explained before. Another possibility would have been to use two power laws which also display two asymptotes. The issue with this option is that there can be a significant overlap of the two power laws that may affect the whole fitting procedure. A two‐terms power law worked well in the case of N2 but due to the scarcer database of O2, it misrepresented the trend that η should have at high energies. With the additional parameter k of a hyperbola (25a), the radius of curvature at the connection point between the two asymptotes can be controlled.

The coefficients in (25) relate to the two asymptotes of the hyperbola through the following: The coefficients thus defined are given for each element of air (N2,O2) in Table 5.

Table 5

Screening Logarithmic Hyperbola Fitting Parameters

	N2		O2
n	a n	b n	a n	b n
1	−2.61	5.24	−1.63	3
2	‐0.54	−1.68	−0.58	−1.31
Murphy	−1.377	5.77	—	—

The major difference comes from the fact that, as mentioned previously, Murphy (1988) used only data from Shyn et al. (1972) (yellow dots on Figure 7, top) to extrapolate his relation. With a more extended database, the power law is not as steep and changes slightly according to each molecule. It appears also that the last yellow point at 90 eV on (Figure 7, top) is an outlier. This could be because the experiments by Shyn et al. (1972) were one of the firsts and only to be able to discriminate vibrational levels of excitation. All data beyond this energy (DuBois & Rudd, 1976; Kambara & Kuchitsu, 1972) include actually inelastic scattering due to vibrational excitation. This is not a great issue as long as the specificities of each experimental data is taken into account. We therefore justify the use of those vibrationally unresolved data by simply applying the same scattering profile for both vibrational and elastic collisions in the Monte‐Carlo model.

Backscattering

A key point to consider when approximating the DCS is the relative contribution of backscattering, which may affect the transport properties of electrons. In Figure 8 we show the DCS against for electrons with 1 keV obtained from Herrmann et al. (1976) as well as the pure‐Rutherford, SR, and modified SR. We clearly see a departure from the Rutherford scattering around 40°. The reason is that the electron‐molecule scattering at intermediate energies and even high energies comprises contributions to the DCS that are not properly accounted by the first Born approximation for screened atomic nuclei. This deviation can be observed at other energies (800, 500 eV) and in the case of O2 as well. It results from a combination of factors, and the following points should be considered.

First, the polarization potential brings a contribution to small‐angle scattering (Christophorou, 2012; Herrmann et al., 1976; Shimamura & Takayanagi, 1984) not included in the SR.

Second, it was already suggested (Hayashi & Kuchitsu, 1976a) that the potential of N2 needs more refinement than that of two decentered homonuclear screened atomic potentials. The chemical bond influence on atomic potentials should be taken into account (Lee & Freitas, 1981). Also, the validity of the first Born approximation, from which the SR is derived, was put into question and more advanced approximations (e.g., eikonal as suggested by  Glauber, 1959) were implemented by Tayal et al. (1983).

Third, intramolecular multiple scattering could affect the DCS shape. The correction provided by including double scattering interference was shown to be negative and that the DCS depended critically on the singly scattered wave function taken (Iga et al., 1999; Kohl & Arvedson, 1980).

Fourth, the hypothesis that the experiment(s) did not properly treat background noise at large‐angle scattering is very unlikely because the results from O2 and N2 both show an underestimation of the DCS by the SR and come from different sources.

In any case, this results in an overall probability shift that affects about 50% of the electrons. The average deviation is overestimated by 5% when neglecting backscattering. This value is not very high on its own, however it introduces a systematic bias: as the electrons experience more collisions, the effective loss of energy will each time be a factor of 1/0.95 higher than it should be. Those errors accumulate with simulation time.

Particle Code

We developed a Monte Carlo code to investigate the effect of DCS for elastic scattering on the transport parameters of an electron swarm in air. The code is fully relativistic and integrates the equation of electrons under given electric and magnetic fields E and B (Jackson, 1998): where p is the kinetic momentum, γ is the Lorentz factor, m is the electron mass, and e is the elementary charge.

The dynamic equation (26) is punctuated by stochastic collisions that change the electron's momentum within a very short time scale that we consider instantaneous. Each electron evolves asynchronously and we consider elastic, excitation, attachment, and impact ionization collisions and implemented a null collision method (Koura, 1986; Lin & Bardsley, 1978) that allows us to use an energy‐independent collision rate at the cost of rejecting some collision events. The total cross sections for all the processes we implemented (including the elastic cross section) come from Phelps and Pitchford (1985) database available on the LXcat server (www.lxcat.net/Phelps). Let us shortly review how we handle each of these collision types.

Elastic Scattering

Whenever an electron suffers an elastic collision under cylindrical symmetry, provided a certain DCS, it is scattered at an angle θ obtained by the cumulative integral: where x is picked from a uniform random variable between 0 and 1, and σ(ε) is the fully integrated DCS (from θ ′=0 to θ ′=π). To invert this relationship in practice, it is possible to store tabulations of the cumulative integrals at the desired scattering precision and retrieve the corresponding angle for each trial. If, however, the DCS has a simple analytical expression (such as the SR (2) or any other, such as the ones in sections 2.2 or 2.3), then (27) can be integrated and inverted to find θ analytically as well. In our case, we use tabulations when including Legendre Polynomials in the DCS but resort to the analytical inverse relation when using SR and other simple formulae.

Excitation

The excitation corresponds to a loss of energy followed by scattering. The cross sections of excitation were taken from Phelps and Pitchford (1985). Technically, excitations have their own differential scattering characteristics. However, when measuring differential scattering cross sections experimentally, it is difficult to discriminate the contribution of inelastic scattering from elastic scattering when the energy resolution of the detector is not fine enough (Bromberg, 1969; DuBois & Rudd, 1976; Iga et al., 1987; Kambara & Kuchitsu, 1972). Therefore, in practice, the experimental DCS for elastic processes can be used as well for most excitation cross sections, which is the approach that we used.

Ionization

The ionization process consists in freeing a (secondary) electron, with a ionization threshold of ε ion, and giving it a certain kinetic energy ε 2. The energy loses sustained by the impacting electron are of the following: where the initial and final energies are denoted ε and ε 1, respectively. Note that since both scattered electrons are indistinguishable, it is customary to assign the highest energy to the primary electron, which implies ε 1>ε 2. We sample ε 2 according to Opal et al. (1971) with the uniformly distributed random variable x∈[0,1]:

Neglecting the kinetic energy transferred to the ionized target, the angles of deviation θ 1 and θ 2 are determined by ε 2 as (Lehtinen et al., 1999)

Results and Discussion

We now apply the code described in the preceding section to investigate the effect of various DCS on electron transport in air, with a focus on the generation of high‐energy electrons. In the text and the figures, we use reduced electric field E/n, which is the electric field divided by the air number density. The runaway threshold, conventional (E ) and critical breakdown reduced fields E/n are independent of the air number density and therefore of the altitude. The traditional unit of the reduced electric field is the Townsend(1 Td = 10−21 V m2) and at atmospheric pressure 1Td corresponds to about 27 kV/m. The conventional breakdown in standard conditions is approximately 3 MV/m = 111.6 Td.

The different scattering schemes are the following:

represents the fit described in section 2.4.

is a linear interpolation of our database extended with the formula by Moss et al. (2006) (12) beyond 1 keV.

is the analytical direct and reversed SR fit (9) provided in Murphy (1988).

is the exponential, direct, and reversed screened fit (3) given by Jackman and Green (1979).

comprises only the SR part beyond 15 eV of this work's fit (16).

is an isotropic distribution in the differential scattering cross section (for comparison).

For each of these choices we ran simulations of electron swarms under homogeneous and constant electric fields ranging from the runaway threshold (∼8 Td or 0.072E ) up to the critical breakdown value (∼1,116  Td or 10E ) (Dwyer et al., 2012). In each simulation the initial swarm comprises 1,000 electrons at an arbitrary temperature of 5000 K. The number of electrons is kept constant during the simulation by randomly dropping one electron whenever an ionization event pushes the number of electrons over this limit. Although it unavoidably increases the noise in the data, this procedure does not introduce a bias because any expected value of a population is the same as that of a randomly sampled subset. In any case, to verify this assertion, numerous comparison tests were performed by comparing big simulations including up to one million particles to those 1,000. No bias could be observed between both batches of simulation. The total simulated time was 1 ns. All macroscopic swarm parameters were observed to converge in less than 0.1 ns at most (lowest electric fields), thereby assuring that equilibrium had been reached. All parameters presented are then averaged in time and space over the remaining 0.9 ns steady‐state part of the simulation. We also compared the results of the Monte Carlo code with two steady‐state Boltzmann solvers: the two‐term Bolsig+ (Hagelaar & Pitchford, 2005) and the multiterm data available from Dujko et al. (2011).

Figure 9 shows the mean electron kinetic energy for different applied electric fields. The range covered starts below 1 eV (runaway threshold) and climbs around 20 eV (critical threshold). We first note that the lines corresponding to linear interpolation (Linterp) and to the fit presented here (SR+LP) practically overlap, which supports the consistency of our constructed DCS (16). The next closest curve is the one using the DCS in Murphy (1988). This DCS was fitted for energies of 5 up to 90 eV, which is why there is very little difference seen in the bulk of the electrons at thermal energies (around their average).

Figure 9

Mean electron kinetic energy as a function of the externally applied electric field for different scattering schemes. Dashed lines are Boltzmann solver codes (as references). See text for the meaning of each label.

Focusing now on electric fields below the conventional breakdown ( 110 Td), we observe that all scattering schemes roughly coincide. The agreement with Bolsig+ and the results by Dujko et al. (2011) (dashed lines in Figure 9) for those fields was used as a validation of our code.

However, for stronger electric fields, the lines diverge. The field enhances, thus the difference in electron scattering profiles implemented at higher energies. The divergence of Bolsig+ calculations from our implementation of the Isotropic scheme is caused by Bolsig+ restriction of the anisotropic velocity distribution in the phase space to an expansion of Legendre polynomial of order 1. This assertion is supported by the better agreement with a multiterm Boltzmann solver (Dujko et al., 2011).

Concerning the comparison between different DCS fits, the one using the DCS in Murphy (1988) intersects and lies close to this work. This DCS was fitted for energies from 5 to 90 eV, which is why there is very little difference seen in the bulk of the electrons at thermal energies. An opposite situation can be seen from the fit provided by Jackman and Green (1979). There is a systematic underestimation of about 10% from the Linterp scheme. As a complement, two extreme DCS cases were provided: an isotropic distribution and a forward scattering distribution represented by only the SR component of the fitted DCS. They provide a background to compare how the DCS in different energy domains affects the swarm spectrum according the the strength of the electric field. At this point, it is crucial to remind that

Jackman‐Green, SR, and Isotropic implementations are all exactly equal for energies below 15 eV. The divergence in the mean electron energy for those three schemes is thus directly related to the dynamics of electrons with energies above 15 eV. This becomes noticeable from around 350 Td which is about 3 times the conventional breakdown field at which the mean energy is 7.5 eV: half the threshold value where the schemes sharply differ. On a macroscopic scale, this is also the point from which it becomes important to model the DCS properly. Neglecting the backscattering term like in the SR case raises the average electron swarm energy, thus its diffusion and drift properties. An isotropic treatment of the scattered electrons reaches its validity limits at fields of about 300 Td. Even a less accurate implementation of the backscattering trend introduces an average bias in the energy spectrum (as in the case of Jackman‐Green).

Finally, it is interesting to note the approximately linear dependency of the average electron energy on the applied electric field in the isotropic case in contrast with the curvatures induced by other schemes. This can be intuitively interpreted as electrons getting proportionally more energy from the field before their next collision occurs. When the scheme is not isotropic for all energies, a change from roughly isotropic to forward scattering regimes causes the electron spectrum to extend farther like in the SR case.

Let us look now at the electron energy distribution function (EEDF) resulting from the different scattering schemes and under a very high electric field. We performed a new set of simulations under the critical homogeneous electric field (1,116 Td) and in this case, starting from a seed of 1,000 particles in steady state taken from the previous simulation outputs, we allowed the particle number to grow up to one million which was set as the upper threshold after which we used the same random filter. The simulation time was at least 0.3 ns and all macroscopic quantities were observed to remain in steady state during the whole simulation. We then averaged the spectra taken over the whole simulation period every 10 ps. Since the field is critical for thermal runaway, somewhat arbitrarily we chose to speak of runaway electrons as electrons with energies beyond 1 keV. The underlying assumption is that the spectral densities at 1 keV are indicative of spectral densities at higher energies in thermal runaway regimes through the balance between the average friction force and the electric force. In a more realistic setting with a highly heterogeneous field around a streamer tip, electrons beyond 1 keV should be considered rather as potential runaway electrons.

The EEDF obtained at equilibrium for the critical field value shows significant differences (Figure 10) according to what DCS is used. We observe a considerable shift in the electron distribution when the SR‐only and Murphy's DCS are used. Their production of runaways is enhanced by more than 2 orders of magnitude. Murphy's curve overruns the SR at 600 eV which can be explained by the much faster drop in its screening parameter than in ours (Figure 7, top). Beyond 1 keV, the LP part in the SR+LP scheme is extinguished (i.e., SR+LP(>1 keV) = SR). However, the distributions yielded between SR and SR+LP do not merge in the graph given with the present resolution. This shows how much the scattering scheme at lower energies affects the EEDF for higher energies. On the other hand, the high‐energy electron production of the Jackman‐Green DCS is not so different from ours as could have been expected; the comparison made in Figure 9 with the average kinetic energy was masking this similarity because a considerable number of electrons is trapped at lower energies in a scheme that is isotropic below 15 eV (in the case of Jackman‐Green). Also, quite unsurprisingly, the use of an isotropic scheme does virtually not allow for runaway production; the spectrum is sharply cut off around 200 eV. Because of the important noise present from 500 eV, it is quite difficult to infer major differences between the runaway output due to Linterp, Jackman‐Green, and our SR+LP curve. In principle, we could observe how a sharp transition in the DCS can stimulate electrons to accelerate further at those energies. For instance, the extension of the Linear Interpolation of the database by Moss' formula accentuates forward scattering which facilitates the production of electrons up to 100 keV in that case. Nevertheless, a more advanced particle management algorithm should be devised to be able to confirm this hypothesis.

Figure 10

Energy spectrum of electron swarms for different scattering schemes at the critical electric field for thermal runaway (1,116 Td). All spectra are normalized (their integral = 1) so that they can be compared. The energy binning is adapted to the logarithmic scale shown on the graph, it is calculated according to Knuth, 's (2019) rule for bins of identical widths in the logarithm of the electron energy domain. The Bolsig+ result (dashed curve) is also provided as an indication although its results do not apply under the critical field.

In the end, two DCS that yield similar bulk properties from their shape around thermal energies ( 50 eV) can lead to very different output in runaway electrons because of their divergence above 100 eV (SR+LP versus Murphy). And conversely, two DCS that can be quite different at low energies can have a similar output of electrons at higher energies (Linterp+Moss and Jackman‐Green).

The relative spectra for each scheme can be qualitatively understood by checking the average deviation committed by a DCS at a certain energy. Its effect can be interpreted as an additional average friction force from elastic collisions. This is represented in Figure 11; for each DCS (the label “x” in the equation below relates to any of the labels in Figure 11), the average was plotted for each energy:

Figure 11

Average scattering‐angle cosine for various DCS fits. The jagged line aspect unveils first‐order discontinuities in the DCS representation discussed earlier. Our work (solid red) strove to smoothen the curve given by pure linear interpolation of the data points (dashed green) starting from 15 eV. A zero mean cosine here implies an isotropic DCS segment.

Here we can observe clearly the weaknesses of some DCS seen as zeroth‐order discontinuities (SR at 15 eV and Jackman‐Green at 15 and 200 eV), and first‐order discontinuities (bends on Murphy's curve). The leaps at 15 eV are due to a sharp switch from an isotropic scheme. At low energy, the most scattering DCS will be able to retain electrons and prevent them from accelerating. This is why, Jackman‐Green's spectrum is the highest below 15 eV (except for the isotropic scheme). Nonetheless, even though the SR is also isotropic beneath 15 eV; the jump in average deviation is so high that it acts as an escape hole: all electrons that by chance reach 15 eV will immediately be accelerated in a runaway fashion. Conversely, in the case of Jackman‐Green, the isotropicity below 15 eV is not masked by the discontinuity and will pull the spectrum toward the isotropic curve until 20 eV when it joins paths with Murphy's curve.

We can also see the extent of the spectrum's bulk influence on the tail: although Murphy is in lead for forward scattering from 300 eV, it is not until 600 eV that its corresponding spectral curve overtakes the SR that had precedence in the number of energetic electrons in the intermediary range. A similar situation happens in comparison with our scheme; it is surpassed in forward scattering by Murphy around 60 to 80 eV and by Jackman‐Green after they cut their backscattering tail at 200 eV, in (5). When we look at the spectra, Murphy's scheme produces more fast electrons from 110 eV. Since the slope of the spectrum governs the balance in electrons being accelerated and decelerated, it is interesting to realize that the slope of Murphy's spectrum starts to differ from ours around 80 eV as to indicate the change in scattering strength. Nonetheless, since Jackman‐Green's scheme scatters more on average than SR+LP until 200 eV, this translates into a lower spectral density of high‐energy electrons. Overall, a sudden change or a persisting difference in the average scattering strength (Figure 11) will have a repercussion on the EEDF (Figure 10) since it alters the balance between the average friction exerted on electrons of neighboring energy categories. In particular, such a difference can be caused by misfitting or simply neglecting the backscattering tail, which can be present even at high energies. This is best illustrated through the graph showing deviation of the DCS from the classical Rutherford scattering (Figure 8).

In this whole procedure we compared different DCS by applying them to identical total elastic cross‐section (as they were reportedly implemented in Monte‐Carlo codes). Notice, however, that cross sections are usually inferred from measurements of the momentum transfer cross section, which depends on the average deviation cosine. This suggests that the total cross sections due to different DCS should all be adapted so as to conserve the momentum transfer cross section rather than the total cross section. This observation is an interesting departing point for a future study.

Conclusion

We have investigated the role of different parameterizations of the elastic scattering DCS for electrons propagating in air, with a focus on swarm parameters and how they affect the generation of runaway electrons. We can summarize our conclusions as follows:

The use of the SR cross section alone underestimates backscattering at all energies and sometimes even forward scattering at very small angles. This underestimation introduces a systematic bias in estimations of the production of runaway electrons.

Above 200 eV, this underestimation is especially sensitive to the energy dependence of the screening parameter η. One case in point is the formula by Murphy (1988) which was derived from only one experimental study (Shyn et al., 1972) with a relatively restricted range of electron energies: 5 to 90 eV. This caused the screening parameter η governing the SR (9) to have a wrong asymptotic trend that strongly underestimated its value with regard to more recent experiments at higher energies (Figure 7). This translates into an overestimation of the forward scattering component, which in turn spurs more runaway production.

On the other hand, the expression by Jackman and Green (1979) was derived based on a more extended database and comprising more detail in the small angle scattering. In this case, however, the SR term had a constant screening parameter (7) acting thus as a tunable background term. Also, although their formula was correctly extrapolated up to 1 keV and beyond, its use is only warranted from 30 eV; for the backscattering term in (8) diverges at 12 eV. Another drawback was that their backscattering term had to be cut abruptly at 200 eV in (5) instead of being gradually damped. Overall, Jackman‐Green's DCS slightly underestimates the scattering angles in the middle‐energy range (several tens of eV) and then compensates by a slight overestimation in the high‐energy range (hundreds of eV). Its validity falls down beyond 10 keV and is thus not appropriate for runaway electron modeling.

The parameterization presented by Moss et al. (2006) implemented beyond 1 keV also overestimates strongly forward scattering (see Figure 6‐bottom), which leads to an overproduction of electron runaways. In their original article, the formula was even used from 500 eV for N2, which must have facilitated the runaway output.

As an alternative, we provide a more accurate formula relying on Legendre polynomials, from 15 eV to 1 keV, that present more flexibility for fitting middle‐ and high‐angle ranges. The extrapolated screening parameter for the SR term was fitted with an extended database from the most recent available experiments.

Finally, it is not clear yet whether an SR‐only extrapolation for DCS is accurate enough in the high energies for modeling the process of thermal runaway. A more reliable analytical expression for elastic scattering that can effectively account for the backscattering tail supported by experimental data might be necessary.

Abbreviations and Symbols

In this section we specify and gather the symbols and names given to various properties and concepts in equations as well as in the text.

ε,v: kinetic energy and velocity vector of an incident electron before being scattered.

θ,φ: respectively, polar and azimuthal angles of scattering (cf. Figure 1).

Ω: solid angle into which the scattered electron is projected (with d   dθ   dφ).

m,−e: electron mass and charge, respectively.

Ze: charge of the target's nucleus.

n: number density of (target) molecules.

η: screening parameter governing forward scattering intensity according to the screened Rutherford expression (2).

σ(ε): total electron‐molecule elastic cross section (function of electron kinetic energy).

: differential electron‐molecule elastic cross section.

A : coefficient applying to the Legendre polynomial of order l in the (polynomial) decomposition of the differential scattering cross section.

Acronyms

Differential (scattering) cross section

Legendre Polynomial (in )

Screened Rutherford (2)

Reversed Screened Rutherford (used in Murphy's fit, cf. (9))

Modified Screened Rutherford (introduced by Moss, cf. (12))

negative EXponential DCS (notation taken from Jackman and Green, cf. (3))

Electron Energy Distribution Function (density of electrons in velocity space)

Linear‐Interpolation (of the DCS)

Table B1

Comparison in Number of Parameters Between Different DCS Analytical Fits; Pseudo‐Parameters as Additional Terms

	Present work	Ivanov	Jackman‐Green	Murphy	Moss
# of parameters	21 + 4	25 + 3	10 + 1	5	1