1 Introduction
Collisions of electrons and ions with neutral atoms are relevant in studies of stellar atmospheres (Mashonkina Reference Mashonkina, Adelman, Kupka and Weiss1996), radio emission in recombination processes of H-II clouds, primordial cosmological recombination of hydrogen (Chluba, Vasil & Dursi Reference Chluba, Vasil and Dursi2010) and plasma fusion (Janev et al. Reference Janev, Langer, Evans and Post1987; Janev, Reiter & Samm Reference Janev, Reiter and Samm2003). Initial studies have focused on scattering and excitation of ground state or low-lying states, in particular for hydrogen atoms. More challenging is the extension to high-lying Rydberg states, for which ab initio quantum calculations become prohibitively untenable. This results in the use of various approximation schemes with unavoidable systematic errors which at times can be, in some observable quantities, of order 100 % (Rolfes et al. Reference Rolfes, Gray, Makarov and MacADAM1993; Nagesha & MacAdam Reference Nagesha and MacAdam2003; Przybilla & Butler Reference Przybilla and Butler2004; Vrinceanu, Onofrio & Sadeghpour Reference Vrinceanu, Onofrio and Sadeghpour2014).
Among the various assumptions in these models, the fact that colliding particles may not share energies according to a Maxwell–Boltzmann (MB) probability distribution, as far as we know, has never been systematically scrutinized in the context of plasma fusion, apart from the possible impact in terms of nuclear fusion reactivities (Onofrio Reference Onofrio2018). This is at variance with the astrophysical plasmas for which deviations from MB have been discussed (Nicholls, Dopita & Sutherland Reference Nicholls, Dopita and Sutherland2012; Nicholls et al. Reference Nicholls, Dopita, Sutherland, Kewley and Palay2013; Storey & Sochi Reference Storey and Sochi2014, Reference Storey and Sochi2015; Draine & Kreisch Reference Draine and Kreisch2018). Deviations from the MB distributions are expected in a tokamak-confined plasma both in the scrape-off layer (SOL), and in the divertor region, due to the lower densities and temperatures experienced by electrons and ions with respect to the core confinement region. The lowest density in these regions 
$(10^{13}{-}10^{14}~\text{cm}^{-3})$ implies that the many-body dynamics is closer to the collisionless regime, and therefore there can be deviations from MB due to lack of thermal equilibration. Moreover, edge localized modes (ELMs) in the SOL region can suddenly release suprathermal particles spoiling a pre-existing MB distribution.
In this work, we undertake a quantitative analysis of the electron-Rydberg hydrogen 
$n$-changing excitation and ionization processes with non-MB distributions, with 
$n$ the principal quantum number. The rate coefficient for generic collisional processes can be written as
$$\begin{eqnarray}k_{if}=\langle v\unicode[STIX]{x1D70E}_{if}\rangle =\int f(v)\;v\hspace{2.22198pt}\unicode[STIX]{x1D70E}_{if}(v)\,\text{d}v=\sqrt{\frac{2}{m}}\int P(E)E^{1/2}\unicode[STIX]{x1D70E}_{if}(E)\,\text{d}E,\end{eqnarray}$$ where 
$\unicode[STIX]{x1D70E}$ is the cross-section for the given process, 
$f(v)$ is the probability distribution of the velocities of the projectile particles and 
$P(E)$ is the corresponding probability energy distribution. The goal is to evaluate changes in the rate coefficients when distributions 
$f(v)$ and 
$P(E)$ in (1.1) deviate from MB.
Deviations from a MB distribution are investigated within two classes, a mixture of two MB distributions at different temperatures, and the so-called 
$\unicode[STIX]{x1D705}$-distribution. These two examples of non-Boltzmann distributions have in common qualitative features, such as the presence of a substantial high-energy tail, and it is therefore interesting to study their quantitative impact with respect to a MB distribution for instance in collisions where the particles share the same low-energy distribution, or total internal energy.
The range of energies and densities is chosen in order to characterize plasmas around the scrape-off layer of a tokamak machine and the divertor region, with electron densities, 
$\unicode[STIX]{x1D70C}_{e}\sim 10^{13}~\text{cm}^{-3}$ and energies, 
$E\sim 0.5{-}20~\text{eV}$ (Anderson et al. Reference Anderson, Ballance, Badnell and Summers2002). We describe in detail the two non-MB distributions in § 2, and then proceed to discuss the results for excitation and ionization of hydrogen atoms of high 
$n$ in §§ 3 and 4, respectively. In the concluding section, we relate the results to the atomic physics in the tokamak divertor region.
2 Non-Maxwellian energy distributions
We have chosen energy distributions differing from MB having in mind examples already available in various physical contexts which possess pronounced, ‘hard’, high-energy tails. This feature can sensibly change the rate coefficients for processes with energy dependent cross-sections. In general, we expect formation of distributions with high-energy tail whenever there is an energy inflow into the system which is large enough to not be dissipated into all the modes during the typical relaxation time scales in the system (Livadiotis & McComas Reference Livadiotis and McComas2011; Livadiotis Reference Livadiotis2018a).
The first example is provided by a mixture of two Maxwell–Boltzmann distributions (bMB) with different inverse temperatures 
$\unicode[STIX]{x1D6FD}_{1}=1/(k_{B}T_{1})$ and 
$\unicode[STIX]{x1D6FD}_{2}=1/(k_{B}T_{2})$
$$\begin{eqnarray}P_{bMB}(E)=2\sqrt{\frac{E}{\unicode[STIX]{x03C0}}}[p_{1}\unicode[STIX]{x1D6FD}_{1}^{3/2}\exp (-\unicode[STIX]{x1D6FD}_{1}E)+p_{2}\unicode[STIX]{x1D6FD}_{2}^{3/2}\exp (-\unicode[STIX]{x1D6FD}_{2}E)],\end{eqnarray}$$ where 
$p_{1}$ and 
$p_{2}$ are the statistical weights of the two distributions, i.e. the probability that a given particle will belong to the distribution with 
$\unicode[STIX]{x1D6FD}_{1}$ or 
$\unicode[STIX]{x1D6FD}_{2}$ (where we assume, for instance, 
$\unicode[STIX]{x1D6FD}_{1}>\unicode[STIX]{x1D6FD}_{2}$), respectively. Since 
$p_{1}+p_{2}=1$, the probability distribution is characterized by three independent parameters, 
$\unicode[STIX]{x1D6FD}_{1},\unicode[STIX]{x1D6FD}_{2}$ and either one of 
$p_{1}$ or 
$p_{2}$. Having in mind cases in which the high-temperature component is not dominant, we choose 
$p_{2}$ such that in the 
$p_{2}=0$ case, the MB distribution is recovered at 
$\unicode[STIX]{x1D6FD}_{1}$.
Such distributions appear in fusion plasma for instance after ion cyclotron resonance heating, as discussed in Bhatnagar et al. (Reference Bhatnagar, Jacquinot, Start and Tubbing1993). The neutral hydrogen flux energy distribution was measured, during heating, to contain a high-energy MB tail with a temperature of 48 keV on top of the pre-existing MB distribution at a temperature of 3 keV. Another example is provided by temperature anisotropy driven instabilities in the solar wind, for which the velocity distribution of the involved particles is adequately fitted by bMB distributions (Klein & Howes Reference Klein and Howes2015; Yoon Reference Yoon2017; Klein et al. Reference Klein, Alterman, Stevens, Vech and Kasper2018). Shocks and winds affect the electron velocity distribution in collisionally ionized plasmas, producing low-energy electrons with power-law tails (Hahn & Savin Reference Hahn and Savin2015; Cui et al. Reference Cui, Foster, Yuasa and Smith2019). The rate coefficients corresponding to these bMB distributions are linear combinations of MB rate coefficients for temperatures 
$T_{1}$ and 
$T_{2}$, with weights given by 
$p_{1}$ and 
$p_{2}$.
More intriguing is the case of the 
$\unicode[STIX]{x1D705}$-distributions. These are generalizations of Lorentzian distributions which were first introduced for applications in space plasma physics (Vasyliunas Reference Vasyliunas1968). The characterization of these 
$\unicode[STIX]{x1D705}$-distributions requires two parameters determining the shape and the corresponding temperature that can be associated with a probability distribution in a sense that is explained below.
The energy probability density is defined by
$$\begin{eqnarray}P_{\unicode[STIX]{x1D705},\unicode[STIX]{x1D702}}(E)=2\sqrt{\frac{E}{\unicode[STIX]{x03C0}}}\unicode[STIX]{x1D6FD}^{3/2}\frac{C(\unicode[STIX]{x1D705},\unicode[STIX]{x1D702})}{\displaystyle \left(1+\frac{\unicode[STIX]{x1D6FD}E}{\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702}}\right)^{\unicode[STIX]{x1D705}+1}},\end{eqnarray}$$ where 
$C(\unicode[STIX]{x1D705},\unicode[STIX]{x1D702})=\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D705}+1)/[\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D705}-1/2)(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702})^{3/2}]$, 
$\unicode[STIX]{x1D705}>3/2$ and 
$\unicode[STIX]{x1D702}>-\unicode[STIX]{x1D705}$. The MB distribution is recovered as 
$\unicode[STIX]{x1D705}\rightarrow \infty$. According to a semiqualitative analysis presented in Livadiotis (Reference Livadiotis2018b), see in particular figure 3, any distribution with 
$\unicode[STIX]{x1D705}>20$ can be well approximated with an equilibrium MB distribution, while typical values for genuine 
$\unicode[STIX]{x1D705}$-distributions out of thermal equilibrium require 
$\unicode[STIX]{x1D705}$ in the range between about 2 and 10. Our choice of 
$\unicode[STIX]{x1D705}$-parameters in the following considerations is based on this criterion.
Notice that we have introduced a parameter 
$\unicode[STIX]{x1D702}$ which allows us to treat various 
$\unicode[STIX]{x1D705}$-distributions differing in their average energy content. Indeed, the 
$\unicode[STIX]{x1D705}$-distribution has a maximum at energy 
$E_{\text{max}}=(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702})/[(1+2\unicode[STIX]{x1D705})\unicode[STIX]{x1D6FD}]$ and an average energy 
$\langle E\rangle =3(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702})/[2(\unicode[STIX]{x1D705}-1)\unicode[STIX]{x1D6FD}]$. This allows for various interpretations of the temperature associated with a 
$\unicode[STIX]{x1D705}$-distribution. For example, by setting 
$\unicode[STIX]{x1D702}=1/2$, the distribution peaks at 
$E_{\text{max}}=1/(2\unicode[STIX]{x1D6FD})=k_{B}T/2$ independently of 
$\unicode[STIX]{x1D705}$ and therefore in practical applications its temperature can be derived as 
$T=2E_{\text{max}}/k_{B}$. For 
$\unicode[STIX]{x1D702}=-3/2$, the average energy 
$\langle E\rangle$ is independent of 
$\unicode[STIX]{x1D705}$, and temperature can be obtained as 
$T=2\langle E\rangle /3k_{B}$. Another interesting case is 
$\unicode[STIX]{x1D702}=0$ for which the 
$\unicode[STIX]{x1D705}$-distribution of velocities has a 
$\unicode[STIX]{x1D705}$-independent maximum corresponding to the most probable velocity 
$v_{\text{max}}=\sqrt{2k_{B}T/m}$ (see appendix A for details). While from an operative standpoint, the values of 
$\unicode[STIX]{x1D702}$ and 
$\unicode[STIX]{x1D6FD}$ can be extracted by looking at the behaviour around the peak of the distribution, the value of 
$\unicode[STIX]{x1D705}$ can be obtained by fitting the high-energy tail.
Absolute (a–c) and relative (d–f) comparison of bMB and 
$\unicode[STIX]{x1D705}$-distributions to MB distributions. The solid lines in the top panels represent the limiting cases for the MB distributions, i.e. 
$p_{2}=0$, and 
$\unicode[STIX]{x1D705}\rightarrow \infty$ for the 
$\unicode[STIX]{x1D705}$-distributions: (a) bMB distributions with increasing weight 
$p_{2}$ of the high-temperature component, and 
$T_{2}=10\;T_{1}$; (b) 
$\unicode[STIX]{x1D705}$-distributions for 
$\unicode[STIX]{x1D702}=1/2$ have maxima at 
$\langle E\rangle =\frac{1}{2}k_{B}T$, indicated by vertical dashed lines; (c) 
$\unicode[STIX]{x1D705}$-distributions for 
$\unicode[STIX]{x1D702}=-3/2$ have average energy given by 
$\langle E\rangle =\frac{3}{2}k_{B}T$, indicated by vertical dashed lines. Panels (d–f) show more clearly the relative deviations from a MB energy distribution for the corresponding top panels, especially in the high-energy tails. Notice that the 
$\unicode[STIX]{x1D705}$-distributions for 
$\unicode[STIX]{x1D702}=1/2$ resemble more closely the bMB distributions, compared to the corresponding 
$\unicode[STIX]{x1D705}$-distributions for 
$\unicode[STIX]{x1D702}=-3/2$. In the latter, the deviations from a MB distribution at high energy are less pronounced, while significant deviations instead occur at low energy. The case of 
$\unicode[STIX]{x1D702}=0$, not depicted for graphical reasons, has a behaviour close to the 
$\unicode[STIX]{x1D702}=1/2$ case.
Figure 1 illustrates the main features of bMB and 
$\unicode[STIX]{x1D705}$-distributions. They all have only one maximum and exhibit strong tails at large energies, particularly manifest in the lower plots, showing the relative difference from the corresponding MB distribution. For 
$\unicode[STIX]{x1D702}=1/2$, the maximum of the distribution is at the same position, showing a significant fraction of particles with higher energy, as compared to MB distributions. At the other extreme, for 
$\unicode[STIX]{x1D702}=-3/2$, the average energy is the same for various values of 
$\unicode[STIX]{x1D705}$, pushing a large number of particles towards energies lower than what would be expected for a MB distribution, thereby counterbalancing the still significant population in the high energy tail.
3 Electron-Rydberg atom excitation rate coefficients
Rydberg atom excitations due to electron scattering have been studied with different techniques and in various physical contexts, ranging from cold atomic plasmas in the laboratory (Rolfes et al. Reference Rolfes, Gray, Makarov and MacADAM1993; Nagesha & MacAdam Reference Nagesha and MacAdam2003) to the primordial cosmological recombination (Chluba et al. Reference Chluba, Vasil and Dursi2010). Since the pioneering experiment of Frank and Hertz, excitation by electron bombardment has been studied extensively for low-lying atomic states. However, there has been much more limited success, in terms of overall accuracy, in the case of Rydberg atomic states. This is not surprising because there is a vast gap between the case of low 
$n$, for which exact quantum mechanical computations are still feasible, and the case of high 
$n$ for which, based on the correspondence principle, classical trajectory Monte Carlo (CTMC) simulations are adequate to describe the processes. More specifically, CTMC calculations (Pohl, Vrinceanu & Sadeghpour Reference Pohl, Vrinceanu and Sadeghpour2008) demonstrated that while previous rate coefficients obtained by Mansbach and Keck (Mansbach & Keck Reference Mansbach and Keck1969) are correct for large energy transfers, significant corrections, singular in 
$1/\unicode[STIX]{x0394}E$, have to be introduced for the proper description of collisions at small energy transfer.
For collisional excitation, the proposed rate formula is (Pohl et al. Reference Pohl, Vrinceanu and Sadeghpour2008)
$$\begin{eqnarray}k_{if}=k_{0}\unicode[STIX]{x1D716}_{f}^{3/2}\left[\frac{22}{(\unicode[STIX]{x1D716}_{i}+0.9)^{7/3}}+\frac{9/2}{\unicode[STIX]{x1D716}_{i}^{5/2}\unicode[STIX]{x0394}\unicode[STIX]{x1D716}^{4/3}}\right]\text{e}^{\unicode[STIX]{x1D716}_{f}-\unicode[STIX]{x1D716}_{i}},\end{eqnarray}$$ where 
$k_{0}=\unicode[STIX]{x1D6FD}\text{e}^{4}/\sqrt{m{\mathcal{R}}}$ (
$k_{0}$ expressed in 
$\text{cm}^{3}~\text{s}^{-1}$ in the cgs system with the electric charge in Gaussian units), 
$\unicode[STIX]{x1D716}_{i}=\unicode[STIX]{x1D6FD}E_{i}$, 
$\unicode[STIX]{x1D716}_{f}=\unicode[STIX]{x1D6FD}E_{f}$, with 
$E_{i}={\mathcal{R}}/{n_{i}}^{2}$, 
$E_{f}={\mathcal{R}}/{n_{f}}^{2}$ the absolute values of the initial and final energies, 
${\mathcal{R}}$ the Rydberg constant and 
$\unicode[STIX]{x0394}\unicode[STIX]{x1D716}=\unicode[STIX]{x1D6FD}(E_{f}-E_{i})$.
Equation (3.1) does not describe correctly the 
$\unicode[STIX]{x1D6FD}\rightarrow 0$ limit, because it has a power-like 
$\unicode[STIX]{x1D6FD}^{s}$ behaviour as opposed to the much slower 
$\log (\unicode[STIX]{x1D6FD})$ dependence expected by the Born approximation. This suggests incorporating the expected Born-like behaviour in the classical formula (3.1) to adequately describe the collision rate coefficients over the whole range of temperatures. By replacing the exponential factor 
$\exp (\unicode[STIX]{x1D716}_{f}-\unicode[STIX]{x1D716}_{i})$ in (3.1) with the ‘quantum factor’ 
$\unicode[STIX]{x0394}\unicode[STIX]{x1D716}\;\unicode[STIX]{x1D6E4}(0,\unicode[STIX]{x0394}\unicode[STIX]{x1D716})$, where we have introduced the incomplete gamma function as
$$\begin{eqnarray}\unicode[STIX]{x1D6E4}(0,\unicode[STIX]{x0394}\unicode[STIX]{x1D716})=\int _{\unicode[STIX]{x0394}\unicode[STIX]{x1D716}}^{+\infty }\frac{\text{e}^{-x}}{x}\,\text{d}x,\end{eqnarray}$$ we obtain an expression for the rate coefficient that has the correct behaviour at both low and high temperatures, and maintains its validity even at large 
$n$. Moreover, the formula may be extended to low 
$n$ by introducing a simple fitting factor that is in the range of unity uniformly across all parameters, and that can be found by direct comparison with the accurate R-matrix results of Pryzbilla and Butler, for transitions between low 
$n$ (Wigner & Eisenbud Reference Wigner and Eisenbud1947; Przybilla & Butler Reference Przybilla and Butler2004).
The resulting expression is (Vrinceanu et al. Reference Vrinceanu, Onofrio and Sadeghpour2014)
$$\begin{eqnarray}k_{if}=k_{0}\left(\frac{\unicode[STIX]{x1D716}_{f}}{\unicode[STIX]{x1D716}_{i}}\right)^{3/2}\left[\frac{22}{(\unicode[STIX]{x1D716}_{i}+0.9)^{7/3}}+\frac{9/2}{\unicode[STIX]{x1D716}_{i}^{5/2}\unicode[STIX]{x0394}\unicode[STIX]{x1D716}^{4/3}}\right]\left(\frac{3.5+0.18n_{f}^{2}}{1+1/\unicode[STIX]{x1D716}_{i}^{5/2}}\right)\unicode[STIX]{x0394}\unicode[STIX]{x1D716}\;\unicode[STIX]{x1D6E4}(0,\unicode[STIX]{x0394}\unicode[STIX]{x1D716}).\end{eqnarray}$$ In Vrinceanu et al. (Reference Vrinceanu, Onofrio and Sadeghpour2014), a comparison between this interpolating formula and many analytical models valid at low 
$n$ is carried out for various transitions and in a temperature range 
$2500<T<250\,000~\text{K}$. Although both plots and tables in Vrinceanu et al. (Reference Vrinceanu, Onofrio and Sadeghpour2014) show disagreement between various models even by a factor two in some cases, it is still worth to explore the effect of bMB and 
$\unicode[STIX]{x1D705}$-distributions within the same model especially considering that the cross-section is quite sensitive to the details of the high-energy population.
The evaluation of the rate coefficients for the bMB distribution is easily expressed as 
$k_{if}^{(bMB)}=p_{1}k_{if}^{(1)}+p_{2}k_{if}^{(2)}$, denoting with 
$k_{if}^{(1)}$ and 
$k_{if}^{(2)}$ the rate coefficients for MB probability distributions at 
$\unicode[STIX]{x1D6FD}_{1}$ and 
$\unicode[STIX]{x1D6FD}_{2}$. In figure 2(a), we plot the rate coefficients for electron–hydrogen scattering versus temperature for the bMB distribution with different weights, all normalized to the analogous rate coefficients for a single MB distribution without the high-temperature component (
$p_{2}=0$).
Rate coefficients for 
$\text{H}(n=10\rightarrow n^{\prime }=11)$ Rydberg excitation as a function of temperature, in the 500–5000 K range, for bMB distribution with different weights for the two components and 
$T_{2}=10T_{1}$ (a) and for 
$\unicode[STIX]{x1D705}$ distributions with various 
$\unicode[STIX]{x1D705}$ parameters and 
$\unicode[STIX]{x1D702}=0$ (b). Rate coefficients are scaled by the corresponding MB rates, obtained for 
$p_{2}=0$ for the bMB distribution, and 
$\unicode[STIX]{x1D705}\rightarrow \infty$ for the 
$\unicode[STIX]{x1D705}$-distribution. Notice the suppression of the rate coefficient for bMB occurring at high temperature with a progressively increasing high-temperature (
$T_{2}$) component, as the high-energy tail of this component becomes ineffective for the transition due to the 
$1/E$ dependence of the collision cross-section.
The evaluation of the rate coefficients for the 
$\unicode[STIX]{x1D705}$-distribution is considerably more involved, requiring extensive CTMC simulations similar to the ones carried out in Pohl et al. (Reference Pohl, Vrinceanu and Sadeghpour2008), but starting from 
$\unicode[STIX]{x1D705}$-distributed configurations. The initial positions and velocities for the Rydberg electrons are generated according to a microcanonical distribution corresponding to the 
$\text{H}(n=10)$ energy level, and a continuous distribution of angular momenta. The initial conditions for the incoming electrons are sampled with random impact parameters and velocities from a 
$\unicode[STIX]{x1D705}$-distribution in (2.2) with 
$\unicode[STIX]{x1D702}=0$. The temperature can then be related to the most probable velocity, independent of 
$\unicode[STIX]{x1D705}$, thereby allowing for a fair comparison between the different 
$\unicode[STIX]{x1D705}$ distributions.
To obtain 
$\unicode[STIX]{x1D705}$-distributed velocities, it is most straightforward to first sample a random number 
$u$ from a Fisher–Snedecor F-distribution (Abramowitz & Stegun Reference Abramowitz and Stegun1972), 
$u\sim F(3,2\unicode[STIX]{x1D705}-1)$, and then assign the velocity as 
$v=v_{T}\sqrt{3u(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702})/(2\unicode[STIX]{x1D705}-1)}$, with 
$v_{T}=\sqrt{2k_{B}T/m}$. Most numerical packages have routines for dealing with the F-distribution (Press et al. Reference Press, Flannery, Teukolsky and Vetterling1992). Alternatively, 
$\unicode[STIX]{x1D705}$-distributed random numbers can be obtained by inverting the cumulative distribution function, as detailed in appendix A.
In figure 2(b), results of the CTMC simulations are also shown, each simulation consisting of 
$4\times 10^{5}$ complete trajectories, with random initial conditions generated according with the rules explained above, classified according to the observed outcome, and repeated for each combination of parameters 
$\unicode[STIX]{x1D705}$ and temperature 
$T$. As 
$\unicode[STIX]{x1D705}$ increases, the probability distribution of velocities for the electrons approaches a MB distribution and the rate coefficients scaled by the corresponding MB values, approach unity.
The most pronounced deviation is observed at low temperatures, when the rate coefficients can be as much as 2.5 times larger than the MB ones. In this specific example, the 
$10\rightarrow 11$ transition for hydrogen corresponds to 28.3 meV, i.e. approximately 280 K. In order to remain in a scattering state, the impinging electron producing the transition should have an energy much larger than the excited electron. This implies that MB-distributed electrons in the low-temperature range around 500 K are disadvantaged with respect to 
$\unicode[STIX]{x1D705}$-distributed electrons with a larger high-energy population. As the temperature is increased, this disadvantage is progressively compensated. In the case of a bMB distribution as in figure 2(a), the use of a high-temperature component progressively results in even lower rate coefficients with respect to a single MB distribution as the temperature is increased. This result is reasonable considering that the collision cross-sections decrease as 
$1/E$, and therefore too large energies for the impinging electrons are ineffective in causing the atomic transitions.
4 Electron-Rydberg atom ionization rate coefficients
We now analyse the case of ionization of the Rydberg atoms, employing the generalized ionization cross-sections for Rydberg atoms introduced in Rost & Pattard (Reference Rost and Pattard1997)
$$\begin{eqnarray}\unicode[STIX]{x1D70E}(E)=(1+1/\unicode[STIX]{x1D6FC})^{\unicode[STIX]{x1D6FC}+1}\frac{E_{M}E^{\unicode[STIX]{x1D6FC}}}{(E+E_{M}/\unicode[STIX]{x1D6FC})^{\unicode[STIX]{x1D6FC}+1}}\unicode[STIX]{x1D70E}_{M},\end{eqnarray}$$ where 
$E$ is the excess energy, i.e. the difference between the absolute energy and the ionization threshold energy, and 
$\unicode[STIX]{x1D6FC}$ is the Wannier threshold exponent (Wannier Reference Wannier1953) characteristic of each target–projectile system at low energies (
$\unicode[STIX]{x1D6FC}=1.127$ for the electron–hydrogen collision). This parameterization of the cross-section has the proper low-energy Wannier threshold behaviour, 
$\unicode[STIX]{x1D70E}(E\rightarrow 0)\sim E^{\unicode[STIX]{x1D6FC}}$, the expected behaviour 
$\unicode[STIX]{x1D70E}\sim 1/E$ for 
$E\rightarrow \infty$, and peaks for 
$E=E_{M}$ such that 
$\unicode[STIX]{x1D70E}(E_{M})=\unicode[STIX]{x1D70E}_{M}$. Notice that the 
$1/E$ asymptotic behaviour is reminiscent of a classical approximation, and therefore does not account for the logarithmic quantum corrections, as would be expected in the Born approximation.
Ionization rate coefficients, normalized to the corresponding MB rates, for various bMB distributions with 
$T_{2}=10\;T_{1}$ (a), and for various 
$\unicode[STIX]{x1D705}$-distributions (b). The dots in the right plot originate from CTMC simulations, including their statistical errors, while the lines are from the numerical integration of (1.1) with the model ionization cross-section of (4.1) taken from the Rost and Pattard (RP) parameterization. The parameters used are: 
$n=10$, 
$\unicode[STIX]{x1D702}=0$, 
$\unicode[STIX]{x1D6FC}=1.127$, and 
$E_{M}=2I_{n}$, where 
$I_{n}$ is the binding energy of the Rydberg atom in state with quantum number 
$n$. The RP rate coefficients for the 
$\unicode[STIX]{x1D705}$ distributions strongly depend on the 
$E_{M}/I_{n}$ ratio since a resonance occurring at higher energies has better overlap with the enhanced high-energy tail of the 
$\unicode[STIX]{x1D705}$ distribution. For instance, by choosing 
$E_{M}=5I_{n}$ we obtain RP rate coefficients larger by a factor 1.8 with respect to the 
$E_{M}=2I_{n}$ case for the 
$\unicode[STIX]{x1D705}=2$ distribution.
In figure 3(a) we plot the ratio between the rate coefficient for a mixture of bMB distributions with different values for 
$p_{2}$, and the corresponding rate coefficient for 
$p_{2}=0$ versus temperature. The calculations were repeated for 
$\unicode[STIX]{x1D705}$-distributions with 
$\unicode[STIX]{x1D702}=0$ and the outcome is depicted in figure 3(b). In addition to the analysis based on Rost & Pattard (Reference Rost and Pattard1997), we have also evaluated the rate coefficients based on the same CTMC simulations performed for the excitation analysis. We have chosen an energy for the peak of the cross-section 
$E_{M}=2I_{n}$ as suggested by numerical simulations (Vrinceanu Reference Vrinceanu2005). This is also in line with various parameterization of experimental data showing that the peak of the ionization cross-section is of the order of 5–10 times the ionization energy for low 
$n$, thereafter approaching a value of 2–3 times the ionization energy at higher 
$n$ quantum numbers (Janev et al. Reference Janev, Langer, Evans and Post1987).
The agreement between the simulations and the analytical interpolation based on Rost & Pattard (Reference Rost and Pattard1997) is remarkably good, and yet not completely surprising because they both share the classical behaviour at high energy, omitting Born-like logarithmic corrections to the cross-section. There are large (up to an order of magnitude at the lowest 
$\unicode[STIX]{x1D705}$) deviations at low temperature, in a range of interest for plasma diagnostics in the scrape-off layer region and in the divertor region of a tokamak.
The effect of a more pronounced high-energy tail becomes minimal at high temperature, as expected for phenomena in which a cross-section with a resonant behaviour appears. However, it is worth remarking that we expect suppression of the rate coefficient in some region of the parameter space for both bMB and 
$\unicode[STIX]{x1D705}$ distributions. This behaviour is emphasized by a side to side comparison for two 
$\unicode[STIX]{x1D705}$ distributions, corresponding to different 
$\unicode[STIX]{x1D702}$ parameters, in figure 4. Notice that for 
$\unicode[STIX]{x1D702}=1/2$, a slight suppression occurs, in the discussed temperature range, only for the 
$\unicode[STIX]{x1D705}=2$ case and at the highest temperature. In the case of 
$\unicode[STIX]{x1D702}=-3/2$ instead the suppression is visible for all values of 
$\unicode[STIX]{x1D705}$, reaching almost a factor of two for 
$\unicode[STIX]{x1D705}=2$ at approximately 2000 K. The suppression is not monotonic, an effect already observed in figure 5 in Nicholls et al. (Reference Nicholls, Dopita and Sutherland2012). This effect is easy to explain in terms of mismatching between the resonance of the cross-section of the process and the energy probability distribution. This is more evident in the case of 
$\unicode[STIX]{x1D702}=-3/2$ and small 
$\unicode[STIX]{x1D705}$, as the energy probability distribution has a large excess at both low and high energies, as visible in figure 1(c,f), i.e. far from the resonance condition for the cross-section. The same suppression is also present in bMB distributions, as barely evident in figure 3(a). In appendix B, we discuss an analytical model based on a convenient parameterization of the resonant cross-section, leading to a simple formula showing enhancement and suppression of the corresponding rate coefficient in the opposite temperature limits.
Evidence for the suppression of the ionization rate coefficients with the 
$\unicode[STIX]{x1D705}$ distributions. The rate coefficients are normalized to the corresponding MB rates, for
$\unicode[STIX]{x1D702}=1/2$ (a) and 
$\unicode[STIX]{x1D702}=-3/2$ (b), for various values of 
$\unicode[STIX]{x1D705}$.
5 Conclusions
We have discussed the sensitivity of the rate coefficients for excitation and ionization in electron–atom collisions to deviations from Maxwell–Boltzmann distributions. The outcome indicates that the rate coefficients may differ significantly from those derived from MB distributions especially at low temperatures, in the range of 1000 K. The results readily translate into rate coefficients for proton-atom collisions, since there is an approximate scaling with the square mass ratio in the excitation case, and analogous parameterization of the cross-section with the Rost and Pattard approach, using different Wannier exponents, for the ionization case.
Our discussion is of interest to tokamak physics for two main reasons. First, an accurate knowledge of the kinetics of the involved ions and atoms in the plasma-edge region is a crucial element for the correct working of ITER and related fusion reactors, and for this reason a dedicated facility, DTT, is under construction in Frascati (Albanese et al. Reference Albanese2017) with the aim of guiding the divertor design for the ITER facility. In the plasma-edge region the plasma is colder than in the core of the reactor, and contains impurities due to the interactions with the materials of the vessel wall. Detailed models of this region are important because these boundary effects regulate the amount of impurities penetrating the core, determining the plasma heat load on the divertor target plates, and therefore the overall performance of the reactor (Reiter Reference Reiter1992; Winter Reference Winter2000). Second, the extraction of the temperature and fluxes for plasma components can be affected by this systematic effect, since so far all calculations for the flux of charge-exchanged fast neutrals escaping from a plasma assume MB distributions (Turaginov et al. Reference Turaginov, Brooks, Buzhinskij, Parker, Pulinetz, Skvortzova and Trofimenko1994; Hollmann et al. Reference Hollmann, Brezinsek, Brooks, Groth, McLean, Pigarov and Rudakov2006).
We believe that the 
$\unicode[STIX]{x1D705}$ distributions may play a major role in the atomic physics around the plasma-edge region, since the conditions for local thermal equilibrium are not easily met due to the steep drop in plasma density. While our study is limited to two specific atomic processes, a more general programme should also be implemented for charge exchange and recombination processes (Takamura et al. Reference Takamura, Ohno, Nishijima and Uesugi2002).
Acknowledgements
D.V. is grateful to Texas Southern University High Performance Computing Center for making the necessary computational resources available, and to the National Science Foundation for the support received through grants PHY-1831977 and HRD-1829184. This work was also partially supported by the National Science Foundation through a grant for the Institute for Theoretical Atomic, Molecular and Optical Physics at Harvard University and the Smithsonian Astrophysical Laboratory.
Appendix A. Cumulative functions for generating 
$\unicode[STIX]{x1D705}$-distributed velocities
Monte Carlo classical trajectory simulations require 
$\unicode[STIX]{x1D705}$-distributed velocities corresponding to the energy probability density (2.2), and we summarize here some of their properties and how to generate them from uniform random numbers. The normalized velocity probability density in three dimensions is defined by
$$\begin{eqnarray}f_{3D}(v;\unicode[STIX]{x1D705},\unicode[STIX]{x1D702})=\frac{4}{\sqrt{\unicode[STIX]{x03C0}}}\left(\frac{\unicode[STIX]{x1D6FD}m/2}{\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702}}\right)^{3/2}\frac{\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D705}+1)}{\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D705}-1/2)}\frac{v^{2}}{\displaystyle \left(1+\frac{\unicode[STIX]{x1D6FD}mv^{2}/2}{\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702}}\right)^{\unicode[STIX]{x1D705}+1}},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D6FD}$ is the inverse temperature and 
$m$ is the particle mass. The parameter 
$\unicode[STIX]{x1D702}$ allows for selecting different classes of 
$\unicode[STIX]{x1D705}$-distributions that have special properties and interpretations, as explained below. The conventional MB distribution is obtained from (A 1) in the 
$\unicode[STIX]{x1D705}\rightarrow \infty$ limit.
The 
$\unicode[STIX]{x1D705}$-distribution in (A 1) peaks at a velocity given by 
$v_{p}=\sqrt{2(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702})/(\unicode[STIX]{x1D705}\unicode[STIX]{x1D6FD}m)}$. This shows that the class of 
$\unicode[STIX]{x1D705}$-distributions with 
$\unicode[STIX]{x1D702}=0$ have the maximum at a
$\unicode[STIX]{x1D705}$-independent velocity. Therefore for this class of distributions the temperature can be defined in a 
$\unicode[STIX]{x1D705}$-independent way as 
$T=mv_{p}^{2}/(2k_{B})$, with 
$v_{p}$ the most probable velocity.
The mean and mean square velocities for 
$\unicode[STIX]{x1D705}$-distributions respectively are
$$\begin{eqnarray}\langle v\rangle =2\sqrt{\frac{2(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702})}{\unicode[STIX]{x03C0}\unicode[STIX]{x1D6FD}m}}\frac{\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D705}-1)}{\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D705}-1/2)}\quad \text{and}\quad \langle v^{2}\rangle =\frac{3(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702})}{2\unicode[STIX]{x1D705}-3}\frac{2}{\unicode[STIX]{x1D6FD}m}.\end{eqnarray}$$ For the class of distributions with 
$\unicode[STIX]{x1D702}=-3/2$ the temperature can be defined in a 
$\unicode[STIX]{x1D705}$-independent way as 
$T=m\langle v^{2}\rangle /(3k_{B})$. For a general 
$\unicode[STIX]{x1D705}$-distribution, the ‘core’ temperature, obtained from the most probable velocity, and the ‘kinetic’ temperature that results from the mean square velocity are different, and depend on the specific 
$\unicode[STIX]{x1D705}$ and 
$\unicode[STIX]{x1D702}$ parameters. Only for the special values 
$\unicode[STIX]{x1D702}=0$, and 
$\unicode[STIX]{x1D702}=-3/2$, does the interpretation of temperature become straightforward. Another important class of distributions is 
$\unicode[STIX]{x1D702}=1/2$ where the most probable energy in the distribution (2.2) relates with the temperature regardless of 
$\unicode[STIX]{x1D705}$, essentially the same as in the MB distribution obtained for 
$\unicode[STIX]{x1D705}\rightarrow \infty$.
The cumulative distribution function related to (A 1) can be calculated exactly as
$$\begin{eqnarray}\displaystyle F(x;\unicode[STIX]{x1D705},\unicode[STIX]{x1D702}) & = & \displaystyle \int _{0}^{x}f_{3D}(t;\unicode[STIX]{x1D705},\unicode[STIX]{x1D702})\,\text{d}t=\frac{1}{x}\sqrt{\frac{\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702}}{\unicode[STIX]{x03C0}}}\frac{\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D705}+1)}{\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D705}+3/2)}\frac{1}{(1+x^{2}/(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702}))^{\unicode[STIX]{x1D705}}}\nonumber\\ \displaystyle & & \displaystyle \times \left[\text{}_{2}F_{1}\left(1,-\unicode[STIX]{x1D705}-\frac{1}{2},\frac{1}{2};-\frac{x^{2}}{\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702}}\right)-1-\frac{2\unicode[STIX]{x1D705}+1}{\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702}}x^{2}\right],\end{eqnarray}$$ where 
$x^{2}=\unicode[STIX]{x1D6FD}mv^{2}/2$, and 
$\text{}_{2}F_{1}$ is the ordinary hypergeometric function. Starting from uniformly distributed random numbers 
$u$ and using the inverse of the cumulative distribution function (A 3) one obtains 
$\unicode[STIX]{x1D705}$-distributed velocities using 
$v=F^{-1}(u)\sqrt{2/(\unicode[STIX]{x1D6FD}m)}$. A practical way to solve the transcendental equation 
$F(x)=u$ is to create a two-way table for 
$x$ and 
$u$, and then use an interpolation to obtain 
$x$ for any 
$0<u<1$.
For one-dimensional problems a 
$\unicode[STIX]{x1D705}$-distribution of velocities has a probability density function
$$\begin{eqnarray}f_{1D}(v;\unicode[STIX]{x1D705},\unicode[STIX]{x1D702})=\sqrt{\frac{\unicode[STIX]{x1D6FD}m}{2\unicode[STIX]{x03C0}(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702})}}\frac{\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D705}+1)}{\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D705}+1/2)}\frac{1}{\displaystyle \left(1+\frac{\unicode[STIX]{x1D6FD}mv^{2}/2}{\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702}}\right)^{\unicode[STIX]{x1D705}+1}}.\end{eqnarray}$$ This distribution has a maximum at 
$v=0$, with zero average 
$\langle v\rangle =0$ and variance 
$\unicode[STIX]{x1D70E}^{2}=\langle v^{2}\rangle =(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702})/(\unicode[STIX]{x1D6FD}m(\unicode[STIX]{x1D705}-1/2))$. The special class of 
$\unicode[STIX]{x1D705}$ distributions with 
$\unicode[STIX]{x1D702}=-1/2$ has a 
$\unicode[STIX]{x1D705}$-independent variance, which leads to a direct interpretation of temperature as 
$T=m\langle v^{2}\rangle /k_{B}$. As expected, for large 
$\unicode[STIX]{x1D705}$ the distribution becomes normal, 
$\lim _{\unicode[STIX]{x1D705}\rightarrow \infty }f_{1D}(v;\unicode[STIX]{x1D705},\unicode[STIX]{x1D702})=\sqrt{\unicode[STIX]{x1D6FD}m/(2\unicode[STIX]{x03C0})}\exp (-\unicode[STIX]{x1D6FD}mv^{2}/2)$. The corresponding cumulative distribution function can be expressed as
$$\begin{eqnarray}\displaystyle F_{1D}(x;\unicode[STIX]{x1D705},\unicode[STIX]{x1D702}) & = & \displaystyle \int _{-\infty }^{x}f_{1D}(t;\unicode[STIX]{x1D705},\unicode[STIX]{x1D702})\,\text{d}t\nonumber\\ \displaystyle & = & \displaystyle \frac{1}{2}+\frac{x}{\sqrt{(\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702})\unicode[STIX]{x03C0}}}\frac{\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D705}+1)}{\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D705}+1/2)}\text{}_{2}F_{1}\left(\frac{1}{2},\unicode[STIX]{x1D705}+1,\frac{3}{2},-\frac{x^{2}}{\unicode[STIX]{x1D705}+\unicode[STIX]{x1D702}}\right),\end{eqnarray}$$ with 
$x^{2}=\unicode[STIX]{x1D6FD}mv^{2}/2$. An alternative method to generate one-dimensional (1-D) 
$\unicode[STIX]{x1D705}$-distributed velocities uses random variates of the Student t distribution, as described in Abdul & Mace (Reference Abdul and Mace2014). This distribution is not useful to generate configurations in three dimensions because the magnitude of the sum of three squares of 1-D 
$\unicode[STIX]{x1D705}$-distributed velocities 
$\sqrt{v_{1}^{2}+v_{2}^{2}+v_{3}^{2}}$ is not 
$\unicode[STIX]{x1D705}$-distributed. This does not happen for MB distributions since a three-dimensional MB distribution is factorizable into three one-dimensional normal distributions.
Appendix B. Enhancement and suppression of scattering rates in non-MB distributions: an analytical example
As visible in figure 2(a) and, to a smaller extent, in the analogous one in figure 3, the rate coefficients for bMB distributions can be enhanced or suppressed depending on the temperature range. Intuitively this depends on the overlap between the energy distribution and the cross-section dependence on energy. In particular, for resonant cross-sections, one expects that at high temperature the overlap in the presence of enhanced high-energy tails is smaller with respect to the corresponding MB distribution thereby suppressing the rate coefficient. In order to show this analytically we discuss the case of a simple cross-section parameterized as
$$\begin{eqnarray}\unicode[STIX]{x1D70E}(E)=\unicode[STIX]{x1D70E}_{M}\left(\frac{E}{E_{M}}\right)^{\unicode[STIX]{x1D706}}\exp [-\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D706}}(E-E_{M})],\end{eqnarray}$$ where 
$E_{M}=\unicode[STIX]{x1D706}/\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D706}}$ is the value of the energy for which the cross-section peaks, and 
$\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D706}},\unicode[STIX]{x1D706}>0$. While this cross-section is not encountered in concrete physical applications, it has the same features of many realistic cross-sections, with a peak value at some intermediate energy and tails falling at lower and higher energies. Its advantage is that it allows us to get a simple expression for the rate coefficient when convoluted with a MB energy distribution at inverse temperature 
$\unicode[STIX]{x1D6FD}_{1}$
$$\begin{eqnarray}P_{MB}(E)=2\sqrt{\frac{E}{\unicode[STIX]{x03C0}}}\unicode[STIX]{x1D6FD}_{1}^{3/2}\exp (-\unicode[STIX]{x1D6FD}_{1}E),\end{eqnarray}$$obtaining a rate coefficient, based on (1.1)
$$\begin{eqnarray}k_{1}=\sqrt{\frac{8}{\unicode[STIX]{x03C0}m}}\frac{\unicode[STIX]{x1D70E}_{M}\exp (\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D706}}E_{M})}{E_{M}^{\unicode[STIX]{x1D706}}}\unicode[STIX]{x1D6FD}_{1}^{3/2}(\unicode[STIX]{x1D6FD}_{1}+\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D706}})^{-\unicode[STIX]{x1D706}-2}\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D706}+2).\end{eqnarray}$$ For a bMB mixture involving another inverse temperature 
$\unicode[STIX]{x1D6FD}_{2}$ and weights 
$p_{1}$ and 
$p_{2}$ we then obtain, by defining the common factor 
$C=\sqrt{8/(\unicode[STIX]{x03C0}m)}\unicode[STIX]{x1D70E}_{M}\exp (\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D706}}E_{M})/E_{M}^{\unicode[STIX]{x1D706}}\unicode[STIX]{x1D6E4}(\unicode[STIX]{x1D706}+2)$
$$\begin{eqnarray}k_{bMB}=C[p_{1}\unicode[STIX]{x1D6FD}_{1}^{3/2}(\unicode[STIX]{x1D6FD}_{1}+\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D706}})^{-\unicode[STIX]{x1D706}-2}+p_{2}\unicode[STIX]{x1D6FD}_{2}^{3/2}(\unicode[STIX]{x1D6FD}_{2}+\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D706}})^{-\unicode[STIX]{x1D706}-2}].\end{eqnarray}$$ By taking the ratio between (B 4) and (B 3), with 
$\unicode[STIX]{x1D6FD}_{2}/\unicode[STIX]{x1D6FD}_{1}=\unicode[STIX]{x1D709}=\text{const.}<1$ the rate coefficient ratio is written as
$$\begin{eqnarray}\frac{k_{bMB}}{k_{MB}}=p_{1}+(1-p_{1})\unicode[STIX]{x1D709}^{3/2}\left(\frac{\unicode[STIX]{x1D709}\unicode[STIX]{x1D6FD}_{1}+\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D706}}}{\unicode[STIX]{x1D6FD}_{1}+\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D706}}}\right)^{-\unicode[STIX]{x1D706}-2}.\end{eqnarray}$$ In the (low-temperature) limit 
$\unicode[STIX]{x1D6FD}_{1}\rightarrow \infty$ we have 
$k_{bMB}/k_{MB}\rightarrow p_{1}+(1-p_{1})\unicode[STIX]{x1D709}^{-\unicode[STIX]{x1D706}-1/2}>1$. In the opposite limit, 
$\unicode[STIX]{x1D6FD}_{1}\rightarrow 0$, the rate coefficient ratio instead tends to 
$k_{bMB}/k_{MB}\rightarrow p_{1}+(1-p_{1})\unicode[STIX]{x1D709}^{3/2}<1$.
