{'fraction_non_alphanumeric': 0.03857965451055662, 'fraction_numerical': 0.01199616122840691, 'mean_word_length': 4.4547199441633225, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bistable solutions for the electron energy distribution function in electron swarms in xenon via Boltzmann equation analysis and particle simulations

At low reduced electric fields the electron energy distribution function in heavy noble gases can take two distinct shapes. This"bistability effect"- in which electron-electron (Coulomb) collisions play an essential role - is analyzed here for Xe with a Boltzmann equation approach and with a first principles particle simulation method. The solution of the Boltzmann equation adopts the usual approximations of (i) searching for the distribution function in the form of two terms ("two-term approximation"), (ii) neglecting the Coulomb part of the collision integral for the anisotropic part of the distribution function, (iii) treating Coulomb collisions as binary events, and (iv) truncating the range of the electron-electron interaction beyond a characteristic distance. The particle-based simulation method avoids these approximations: the many-body interactions within the electron gas with a true (un- truncated) Coulomb potential are described by a Molecular Dynamics algorithm, while the collisions between electrons and the background gas atoms are treated with Monte Carlo simulation. We find a good general agreement between the results of the two techniques, which confirms, to a certain extent, the approximations used in the solution of the Boltzmann equation. The differences observed between the results are believed to originate from these approximations and from the presence of statistical noise in the particle simulations.


Introduction
The physics of electron swarms has been attracting considerable attention during the past decades because of the interest in the various physical effects taking place in these settings and by the need for accurate input data in discharge modeling [1,2]. The advance of the theoretical background, as well as of the computational resources and numerical techniques made it possible to develop a detailed picture of the physics of particle swarms. Most of the efforts have been devoted to electron swarms, the electron energy distribution function and the transport properties have been determined for a wide variety of gases and gas mixtures, and for a broad range of conditions, see e.g. [3][4][5][6][7][8][9][10][11][12][13]. One of the particular phenomena, the bistability of the electron energy distribution function (EEDF) is a pronounced manifestation of the nonlinear nature of these systems. The term "EEDF bistability" is used here to designate a situation when in a physical system, under fixed conditions (gas number density, gas temperature, electric field strength, electron and ion concentrations, population of excited states), two stable steady states are possible with different EEDFs. The electron temperature bistability was (according to our best knowledge) predicted for the first time for a plasma under the effect of an alternating electromagnetic field, at conditions when electron-ion collisions are dominant [14,15]. Assuming a Maxwellian EEDF, the bistability in the plasma of heavy rare gases (argon, krypton and xenon) with an applied electric field, was shown to exist in [16] due to the specific shapes of the momentum transfer cross sections for the electron scattering from Ar, Kr and Xe atoms. Later, the effect in heavy rare gases plasma was confirmed by using a more accurate approach, the Boltzmann equation (BE) analysis that took into account electron-electron (e-e) collisions [17][18][19]. In [17][18][19] it was found that, within a certain range of parameters: (i) the reduced electric field strength, E/n and (ii) the electron to neutral atom density, n e /n (where E is the electric field strength, n is the atom number density, and n e is the electron density), the BE has two stable solutions. In the case of xenon gas, e.g., the bistability effect was found at n e /n ≥ 10 −9 and 0.025 Td < E/n < 0.043 Td. The conditions mentioned above can be realized in different physical systems, e.g., in decaying plasmas in the presence of low electric field, in non-self-sustained discharges, and under electron swarm conditions. A decaying plasma in Xe was studied in Ref. [16], where the time-dependence of the current in the afterglow plasma in the presence of a weak electric field was measured. A jump-like decrease in the current was observed at some time point during the plasma decay, which was attributed to the manifestation of the bistability effect [16]. Non-self-sustained discharges (sustained by a beam of fast electrons) were considered in Refs. [20,21], in which the possibility of the existence of the EEDF bistability was theoretically analyzed in Ar, Kr and Xe. It was shown that the bistability effect takes place in Xe and Kr, while for an Ar plasma [21], the BE has a unique solution over the entire parameter range examined. To our best knowledge, the bistability effect has not been investigated so far for electron swarm conditions. We note that besides pure noble gases, the bistability effect has also been studied in gas mixtures. In [22] the electron temperature (T e ) in Ar:N 2 = 100:1 mixture afterglow was studied both experimentally and theoretically. In the experiments it was observed that, under certain conditions, a rather sharp knee appears in T e (t), while the vibrational temperature of nitrogen molecules, T v , remains almost constant. The theoretical study of the EEDF in Ar:N 2 afterglow plasma was carried out by the numerical solution of the appropriate BE, by taking into account e-e collisions, as well as superelastic vibrational and superelastic electronic collisions. The calculations have shown that ranges of n e and T v exist, where two different solutions of BE can be obtained. The observed knee-like time-dependence of T e was explained as the manifestation of the bistability effect. The possibility of this effect in pure nitrogen afterglow plasma was theoretically addressed in [23]. In addition, using the Maxwellian distribution function approach the bistability effect was studied in [24] for positron swarms in He. For more details the reader is referred to the review [25]. It should be noted that so far the theoretical analysis of the bistability effect was performed based on BE approach only. In these studies, and generally, in the solutions of the Boltzmann equation including electron-electron collisions, rather serious approximations have traditionally been used: the methods (i) search for the distribution function in the form of two terms ("two-term approximation"), (ii) neglect the electronelectron part of the collision integral for the anisotropic part of the distribution function, (iii) treat Coulomb collisions as binary events, and (iv) truncate the range of the electronelectron interaction beyond a characteristic distance. These approximations became widely accepted throughout the years, and one needs to note that in the absence of a more rigorous approach, their effects on the results of BE solutions have never been critically examined. To do that -and this is actually the motivation of our work -one would ideally describe the physical system of interest with a first-principles method that does not involve any approximations. (We note that previous particle-based (Monte Carlo) methods, dealing with Coulomb collisions, are neither free of approximations.) The particle-based simulation technique, used here [26], on the other hand, relaxes both the above major approximations and provides a solution based on first principles. Therefore we expect this method to be able to evaluate the accuracy of the BE solution that is limited by the approximations adopted, and can verify the effects predicted within the frame of the BE approach, like the EEDF bistability that was introduced above and is studied here in details in Xe gas. In section 2 we introduce the underlying physical processes responsible for the development of the bistability effect. This is followed in section 3 by the presentation of the methods of calculation: section 3.1 discusses the particle simulation method, while section 3.2 outlines the solution of the Boltzmann equation. The results obtained by the two different approaches are presented and compared to each other in section 4. A short summary is given in section 5.

Qualitative explanation of the EEDF bistability effect in rare gases
The qualitative explanation of the possibility of the EEDF bistability effect in rare gases at weak electric fields (where only elastic collisions occur) adopts an approximation of a Maxwellian EEDF (i.e., it assumes a priori that the EEDF is Maxwellian, f 0 (ε) ∼ exp(−ε/k B T e ), for details see e.g. [25]). For such conditions, the time-dependent equation for the electron temperature is as follows: where the term H(T e ) describes the heating of electrons by the electric field and the term L(T e ) accounts for the loss of electron energy in elastic collisions. (Note that, generally, the loss rate is positive, as electrons deposit energy in collisions with the background gas. However, at T e < T g , where T g is the gas temperature, electrons are heated by the atoms of the gas, resulting in a negative loss, L < 0.) The number of steady state and stable solutions of (1) depends on the shape of the Φ(T e ) function. Figure 1(a) shows this function together with the H(T e ) and L(T e ) functions [25], for Xe at E/n = 0.035 Td, and a gas temperature of T g = 300 K. Let us note that L(T e ) is proportional to the elastic momentum transfer cross section σ m (ε) of electron scattering from Xe atoms, while H(T e ) varies proportionally to 1/σ m (ε) (see section 3a). Due to the specific shape of σ m (ε) (a deep Ramsauer-Townsend minimum at an energy of ≈ 0.6 eV, see figure 1(b)) the H(T e ) and L(T e ) curves exhibit a "wavy" behavior and cross each other at three points. As a consequence, the Φ(T e ) function has an "inverse-N-type" shape and is equal to zero at three T e values. Among these three steady state solutions of equation (1) two (where Φ(T e ) has a negative slope) are stable, as indicated in figure 1(a). Although the EEDF in reality is generally not Maxwellian, but the e-e collisions drive the EEDF towards Maxwellian, therefore the above explanation holds at least qualitatively. We note that e-e collisions are essential in establishing the bistability effect, as the terms corresponding to these collisions introduce non-linearity into the BE. Without e-e collisions the BE is linear equation (in f 0 ). As a linear equation, it has a unique solution at fixed parameters (taking into the account the proper normalization and the condition that f 0 → 0 at electron energy → ∞). If e-e collisions are involved, the BE becomes a non-linear equation and, in principle, two or more different solutions are possible. For multiple solutions to appear the e-e collision frequency should be high enough to have an influence on the shape of the EEDF, i.e., it should be comparable or higher than the frequency of energy losses in electron-atom elastic collisions. We note, however, that the non-linearity is a necessary, but not a sufficient condition for BE to have multiple solutions. While some examples of physical systems, where bistability is present, were given in section 1, a general sufficient condition for the appearance of the effect cannot be formulated at present.

Methods
Here we describe the basics of the two computational approaches, the particle simulation scheme is presented in section 3.1, while the solution of the Boltzmann equation is discussed in section 3.2.

Particle simulation
A swarm of electrons is simulated under the effect of a homogeneous electric field, in xenon gas, at different values of the number density ratio, η = n e /n. At the low reduced fields considered only elastic collisions take place. As this ensures conservation of the number of particles, we consider a spatially homogeneous system. The simulation scheme is based on a combination of a Molecular Dynamics (MD) technique and a Monte Carlo (MC) approach [26]. The MD describes the many-body interactions (driven by the inter-particle Coulomb potential) within the classical electron gas, while the MC part handles the interaction of the electron gas with the background (atomic) gas.
3.1.1. Molecular Dynamics simulation of the electron gas. In our classical many-body system the particles move under the influence of interparticle and external forces in a cubic simulation box, with periodic boundary conditions. The equations of motion of particles i = 1, . . . , N, which are numerically integrated with discrete time steps (∆t), are Newtonian: where m and −e are the mass and the charge of the electron, respectively. The sum on the right hand side represents the force exerted on particle i by all other (j = i) particles and their periodic images located in spatial replicas of the simulation box. These images have to be included in the proper determination of the interparticle forces in the case of the "full" (un-truncated) infinite-range Coulomb potential. (Note that for our conditions no screening of the potential takes place by oppositely charged species.) This summation is indeed a key issue and needs a special approach, as will be explained below. The second term on the right hand side of the above equation is a contribution due to the external electric field. We note that in the absence of the interaction of the electrons with a background gas the electrons would continuously be accelerated by this field, which will, however, not be the case when e − +Xe atom collisions take place.
The MD method adopted in this study is based on the Particle-Particle Particle-Mesh (PPPM) approach, described in details in [28], that uses a partitioning of the interaction into (i) a force component that can be calculated on a mesh (the "mesh force") and (ii) a short-range ("correction") force, which is to be applied to closely separated pairs of particles only. In the mesh part of the calculation, charge clouds are used instead of point-like charges. The charge density distribution ρ(r) is assigned to a grid and is Fourier transformed to the k-space. Multiplying ρ(k) with an optimized Green function results in a potential distribution φ(k) = G(k)ρ(k), which is subsequently transformed back to real space. The forces acting on the particles are obtained by differentiation of the potential and interpolating the electric field to the positions of the particles. The cloud shape is chosen in a way to ensure that ρ(k) is band-limited (which would not hold for point-like charges). The calculation of the potential in the Fourier space automatically takes into account the periodic images of the primary computational box. For closely separated particles a correction force is to be applied, which is the difference between the forces between two point-like charges and two charges with the cloud shape used. The simulations describe a micro-canonical ensemble, where the number of particles, the volume of the system and energy are conserved. The upper limit for the simulation time step is defined by the stability of the integration of the equation of motion in the case of the closest approach of two electrons, r min = e 2 /(4πǫ 0 ε max ). Here ε max is a pre-defined maximum energy [28], which has to be chosen carefully, to ensure that the probability of finding electrons with ε > ε max is vanishingly small at the conditions considered. In this work we adopt a maximum energy of 5 eV, for this value the time step has to be chosen to be as low as ∆t ∼ = 2.2 × 10 −16 s. This results in a very demanding computational load to follow the evolution of the system for a long time. As a converged solution assumes sufficient interaction between the electrons and the gas as well, one needs to have a high-enough electron-atom collision frequency. This can only be ensured by setting a high the gas atom number density.
For this reason our computations are carried out with n = 2.5 × 10 27 m −3 , which allows simulated times of tens of nanoseconds (with a computational speed of simulating ∼ 1 ns in one day, on a single CPU). The number of electrons is chosen to be N = 1000.
3.1.2. Monte Carlo simulation of electron gas -background gas interaction. Having solved the description of the electron gas with the MD method that accounts for electronelectron interactions, now we introduce a background gas and let the two gases to interact via e − +Xe collisions. The probability of an e − +Xe collision during a time step ∆t is calculated as: where g = |g|, with g = v − V being the relative velocity between the electron and a Xe atom with a velocity V randomly chosen from a Maxwellian background of gas atoms having a temperature T g . This probability is calculated for each electron in each time step, and decision about the occurrence of a collision is made by comparing P coll with a random number. Collisions are executed in the center-of-mass frame, and are considered to be isotropic. The energy change of the gas atoms colliding with the electrons is not accounted for, the gas temperature is kept constant. The simulations allow investigation of both the stationary state and the transients (induced by changing the electric field strength or the gas temperature). The transients are followed by monitoring the time-dependence of the mean energy of the electrons, which is calculated in each time step. The EEDF is, on the other hand, only calculated for the stationary state as its accurate determination (over 6-7 orders of magnitude) requires averaging over millions of time steps.

Solution of the Boltzmann equation
Here we discuss the specific features of the BE for electrons, as applied to conditions under consideration (homogeneous plasma, or an electron swarm in an atomic gas, acted upon by a weak steady electric field). The distribution function of the electrons, f (v), can be described by the equation where v is the electron velocity and C is the collision integral. Further we shall consider only elastic scattering of electrons from atoms and electron-electron collisions: The solution of the BE is based on expansion of the distribution function in Legendre polynomials P n (cos Θ), in which only two first terms are taken into account: where v is the velocity magnitude, Θ is the angle between v and −E, f 0 (v) is the symmetrical part of the distribution function and f 1 (v) describes the directed motion of the electrons along the electric field. (We note that this "two-term approximation", in the absence of e-e collisions, has been benchmarked with other solution methods of the BE, as well as with particle-based (Monte Carlo) simulations in several studies, see e.g. [29].) The substitution of expansion (5) into equation (4) leads to equations for the f 0 and f 1 functions: and The collision integrals C 0m and C 1m can be written as [15]: where ν m = nσ m v is the momentum transfer frequency and δ = 2m/M is the average fraction of the energy lost by the electrons in one elastic collision with atom (M is the mass of the gas atom). The rate of the electron energy loss due to elastic collisions is characterized by the frequency ν u = δν m . It is known that the calculation of the pair-collision frequency in the case of Coulomb collisions encounters a characteristic difficulty, namely the logarithmic divergence of frequency at small scattering angles. This difficulty is avoided by assuming that the Coulomb potential acts only up to a certain finite distance r max (see later). The expression for the term S 0e is written as follows [15,30]: where f 0 (v) is normalized as: and ν e = 2πn e e 2 4πε 0 m In BE calculations the value of the parameter r 0 is usually estimated by means of the mean electron energy, that gives r 0 = e 2 /4πε 0 2ε . For the case of plasmas the Debye length, λ D , is taken to be the cutoff distance, i.e. r max = λ D . For the case of swarm conditions considered here we use the approximation: i.e. the half of the average distance between the electrons is taken to be the cutoff distance. (This choice of r max is based on intuitive physical considerations: if the impact parameter (of test electron relative to a given electron) is higher than the half the average distance between electrons in the gas, then the influence of this electron on the test one becomes weaker than the influence of the neighboring one.) Note that in most cases r max /r 0 ≫ 1 (see comments in [15,30]) and the "1" in the expression under the logarithm in (14) can be omitted. As to the term S 1e in equation (4), it is very complex (see comments in [15]) and we did not find publications in which this term was taken into account in calculations. As a rule (BOLSIG+ [31], EEDF [32,33]), it is neglected assuming that ν e ≪ ν m . Then, if the characteristic time of plasma parameters variation is essentially smaller than ν −1 m , the time derivative in equation (7) can be omitted. In this case and equation for the f 0 function is written as It should be noted that, in the absence of e-e collisions, the parameter for the steady state solution of equation (17) is the reduced electric field E/n. If the e-e collisions are taken into account there are three parameters: E/n, n e /n and n e . The electron number density is an independent parameter since the logarithmic term in eq. (14) (the Coulomb logarithm) depends on n e . Actually, at fixed E/n and n e /n values the dependence of f 0 on n e is rather weak, since n e value is under logarithm.
For the numerical solution of (16) it is rewritten with energy as a variable. The steadystate equation is solved by an iteration method similar to that in [33,34]. The initial f 0 (ε) is assumed to be the Maxwellian with a given electron temperature T e0 . In the case of calculation of the time-dependent solution of the BE the time step should be as small as ∆t ≪ ν −1 u and ∆t ≪ ν −1 e . For conditions under consideration ∆t was taken as ∆t = 2 × 10 −14 s.

Results
Bistability of the EEDF, i.e. two solutions for the distribution function, for certain conditions are found here, both via the solution of the Boltzmann equation and via executing the particle-based simulations. To find these different solutions both methods start from different initial conditions: both assume a Maxwellian distribution with a "low" and a "high" mean energy. For the particle methodε init = 0.009 eV and 0.9 eV are used as "low-energy" and "high-energy" initial conditions, with random initial positions of the particles in the simulation box. For the BE solution the initial mean energy values are:ε init = 0.065 eV and 0.65 eV. In the domain where a unique solution exists the results of the calculations do not depend on these choices. In the case of two possible solutions two domains ofε init exist: from one of these domains (ε init <ε * init ) the calculations converge to the lower energy solution, and from the other domain (ε init >ε * init ) to the higher energy solution. The value ofε init within these domains has no effect on the results. The boundary value,ε * init , depends on the E/n and η. The BE analysis has been carried out for a range of reduced electric fields, 0.015 Td ≤ E/n ≤ 0.05 Td, for various electron to gas density ratios. Results of these calculations, in terms of the mean electron energyε, for η = 10 −8 , 10 −7 , and 10 −6 are displayed in figure 2(a). Two solutions are found for the intermediate values of E/n, the boundaries of the domain change slightly with the density ratio. The "low-energy" solution exhibits a mean energy that is almost independent of η, and is in the rangeε ∼ 0.04 . . . 0.05 eV. In the case of the "high-energy" solutions the mean energy amounts several tenth of an electron Volt and increases with increasing η. The particle simulation results agree generally well with those obtained from the BE, small differences are found in the case of "high-energy" solutions. Additionally, the particle simulation does not predict a bistability at η = 10 −8 at E/n= 0.035 Td, whereas the BE solution does. The reason of this discrepancy is not fully understood. We contemplate that the presence of noise in the particle simulation, due to its statistical nature, prevents finding solutions where the minimum in the energy balance vs. T e is very shallow. This, however, needs further clarification. While all the data in figure 2(a) have been shown for a gas temperature of T g = 300 K, figure 2(c) shows the dependence of the results on T g , for fixed η = 10 −7 . The mean energy here is displayed on a log scale, as the gas temperature influences predominantly the "low-energy" solutions. For a lower T g a lowerε is found, however, the domain of E/n, where bistability is found, is wider at T g = 0 K. The stationary solutions in the particle simulations were typically obtained beyond a few tens of nanoseconds of simulated time. The convergence of the mean energy is illustrated in figures 3(a)-(c), for reduced electric fields of E/n = 0.025 Td, 0.035 Td, and 0.045 Td, respectively, at η = 10 −7 . In all cases the simulations were started from two initial configurations, as already mentioned above. While at E/n = 0.025 Td and 0.045 Td we observe convergence to a unique value ofε, the runs at E/n = 0.035 Td clearly yield two solutions, with mean electron energies differing by almost an order of magnitude. The full EEDFs, obtained by the two methods, are compared in figures 3(d)-(f). The EEDFs are normalized as ∞ 0 f (ε) √ εdε = 1. We find a good agreement between the data obtained from the BE and from the particle simulation, although some differences show up in the tails of the distribution functions. The difference of the EEDFs belonging to the two stable solutions at E/n = 0.035 Td is remarkable. To demonstrate the ability of our methods to follow the temporal evolution of the system here we take as an example the change of the electric field. We start from the two stationary stable solutions obtained at E/n = 0.035 Td (at η = 10 −6 ), change E/n to a "new" value, and monitor the convergence ofε(t). We test "new" values of E/n = 0.025 Td and 0.045 Td. The results obtained by both methods are depicted in 4. Following the change of the reduced field both methods -although with a different dynamics -show the relaxation of the system to unique states, corresponding to the new values of E/n. The relaxation time is found to be of the order of few times ten nanoseconds. Here we recall the importance of the energy loss frequency ν u , as this defines the time scale of the relaxation of the mean energy. An estimation at ε = 0.5 eV gives ν −1 u ≈ 20 ns, in accordance with our observations of the relaxation time scale. The peculiarities of the transitions can be understood by following the heating and the cooling of the electrons at the specific conditions. The rate of electron heating is proportional to E 2 /σ m (ε) and the rate of energy loss is proportional to σ u (ε) ε = (2m/M)σ m (ε) ε. First we discuss the transition from the "low-energy" solution at E/n = 0.035 Td to the final state at E/n = 0.045 Td. The change of the mean energy in this case is slow at the beginning but becomes very fast afterwards. In the initial state the EEDF is narrow and the mean energy is low, where the momentum transfer cross section is high. As the heating rate is proportional to σ −1 m (ε) and the cooling rate is proportional to σ m (ε), the rate of change ofε is low at high σ m in the first phase of the transition. When the mean energy becomes higher σ m (ε) decreases, giving more preference to the heating. This, together with the increasing slope of σ m (ε), results in an abrupt increase ofε beyond a certain time.
The transition from the "high-energy" solution at E/n = 0.035 Td to the final state at E/n = 0.025 Td is more complex. The three stages of relaxation -seen in figure 4 -can be explained as follows. Under the steady state conditions at E/n = 0.035 Td the mean electron energy is relatively high,ε ≈ 0.48 eV. At this energy the momentum transfer cross section is low (due to the Ramsauer-Townsend minimum) and, as a consequence, the rate of electron heating is relatively high. The high rate of heating is balanced by a high rate of energy loss. The decrease of the electric field causes an instant decrease of the heating. The unbalanced high rate of energy loss leads to the rapid decrease in the mean electron energy (during the first stage of ε(t) relaxation) and, consequently, to a decrease of the rate of energy loss. The decrease in the rate of energy loss leads to near-balanced conditions, i.e. the difference between the rate of energy loss and the rate of heating becomes moderate, leading to the flattening of the ε(t) curve (second stage of relaxation). When the mean energy decreases below ε ≈ 0.2 eV, the momentum transfer cross section increases sharply and the rate of energy loss increases noticeably, while the rate of heating decreases significantly. Consequently, ε decreases faster below 0.2 eV, during the third stage of the relaxation. We note that since under the conditions considered here the electron-electron collision frequency is higher than the electron-atom energy exchange frequency (ν u ), the EEDF is nearly Maxwellian during the relaxation process. The transitions, which involve a small change of the mean energy (from the "highenergy" solution at E/n = 0.035 Td to the final state at E/n = 0.045 Td, and from the "low-energy" solution at E/n = 0.035 Td to the final state at E/n = 0.025 Td) are much faster, compered to the two cases discussed above.

Summary
We have investigated the bistability of the EEDF in Xe gas at low reduced electric fields via the solution of the Boltzmann equation and via a first principles particle simulation technique. The solution of the Boltzmann equation adopted the usual, widely accepted approximations: it (i) searched for the distribution function in the form of two terms, (ii) neglected the electron-electron part of the collision integral for the anisotropic part of the distribution function, (iii) treated Coulomb collisions as binary events, and (iv) truncated the range of the electron-electron interaction beyond a characteristic distance. The particle simulation method [26], being devoid of any of these approximation has provided first-principles solutions to the problem, via a combination of a Molecular Dynamics simulation method (that described accurately the many-body interactions within the electron gas governed by the full Coulomb potential) and a Monte Carlo method (that handled the interaction of the electrons with the atoms of the background gas). Both methods allowed the computation of the EEDF and the related quantities, and have indicated the existence of two stable solutions for the EEDF for a range of E/n. The electron mean energies and the full EEDFs, obtained by the two methods agreed generally well for most od the parameter settings covered. Differences found for the domain of the bistability and for the shapes of the EEDFs may be attributed to approximations adopted in the BE solutions and to the presence of statistical noise in the particle simulations.